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1.   INTRODUCTION 

In  a  large  supermarket  retail  chain  with  several  storing 
locations,  one  of  the  problems  in  supplying  the  individual 
stores  is  transportation,  i.e.,  how  much  of  a  particular 
item  should  be  sent  each  time  and  from  what  location  should 
it  be  sent  in  order  to  satisfy  the  demands  of  all  the  stores, 
given  the  resources  available  at  different  storing  locations. 

Optimization  of  a  transportation  problem  with  linear 
cost  function  can  be  regarded  as  a  generalization  of  the 
assignment  problem,  which  can  be  carried  out  by  means  of  the 
simplex  method  of  linear  programming  (4).   However,  some 
special  methods,  such  as  the  northwest  corner  method,  the 
unit  penalty  method  and  Vogel's  approximation  method,  have 
been  developed  which  are  easier  to  apply  and  are  less  tedious 
than  the  simplex  method  (1,  11).   Optimization  of  transpor- 
tation problems  with  non-linear  cost  functions  cannot  be 
solved  by  linear  programming  methods.   The  non-linear'  problem 
with  two  and  three  depots  has  been  solved  by  dynamic  program- 
ming (2).   Recently,  a  discrete  version  of  the  maximum  prin- 
ciple has  been  applied  to  the  two-depot  problem  (5,  6)  and 
the  three-depot  problem  (7).   This  has  resulted  in  a  great 
simplification  in  numerical  calculations. 

The  purpose  of  this  report  is  to  study  and  understand 
more  thoroughly  the  application  of  the  discrete  version  of 
the  maximum  principle  in  the  transportation  problem  with  non- 
linear cost  function  so  as  to  enable  us  to  modify  the 


algorithm  in  Ref.  (7)  and  to  write  computer  program  according 
to  this  algorithm. 

In  section  2  the  transportation  problem  is  defined  and 
a  table  of  transportation  costs  and  requirements  is  pre- 
sented.  The  modification  of  the  discrete  version  of  the 
maximum  principle  is  developed  to  derive  the  Hamiltonian 
function  for  the  transportation  problem  in  section  3-   Also 
in  this  section  the  numerical  method  in  obtaining  the  op- 
timal solution  is  illustrated.   Section  4  deals  with  the 
explanation  of  the  steps  in  the  computer  program. 


2.      STATEMENT  OF  THE  PROBLEM 

Suppose  that  there  are  a  number  of  origins  (depots) 
where  the  resource  is  located  and  N  number  of  destinations 
(demand  points)  where  the  demand  for  this  resource  exists. 
Also  suppose  that  there  is  only  one  type  of  resource  and 
that  its  total  supply  is  equal  to  the  total  demand.   Let 

6n  =  the  quantity  of  the  resource  sent  from  the 

i-th  depot  (origin)  to  the  n   demand  point,  and 
Fn(61?)  =  the  cost  incurred  by  this  operation. 

If  there  are  s  depots  and  N  demand  points,  the  problem  is  to 
determine  the  values  of  8. ,  i  =  1,  2,  ...,  s;  n=l,  2,  ..., 
N,  so  as  to  minimize  the  total  cost  of  transporting  the 
resources, 

csN-  I    t    Fj(eJ)  (i) 

n=l  i=l 

subject  to  the  constraints. 

(i)  e">o, 


(ii)  L  6-  =  ™j i  number  of  units  of  the  resource  • 
n=l  x    x 

available  at  the  i-th  depot,  i=l, 

2,  ...,  s, 


(iii)   L.     6  =  D  ,  number  of  units  of  the  resource 

required  by  the  n   demand  point, 

n=l,  2 N. 

A  nonlinear  cost  function  of  the  transportation  problem 
is  as  follows: 

F?(6i)  =  ai9i  '    bi(9i)2  <2)* 

v;here  a.,  b.  are  constants.   The  values  of  a1},   b1}   with  Dn, 
11  l '   l        ' 

the  number  of  units  of  the  resource  required  by  the  n^ 
demand  point,  and  W. ,  the  number  of  units  of  the  resource 
available  at  the  ith  depot,  are  shown  in  Table  1. 


*The  superscript,  n,  indicates  the  stage  number.  The 
exponents  are  written  with  parentheses  or  brackets  such  as 
(e")2  or  t<f>(en)}2. 


Table  1.   Transportation  costs  and  requirements 
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3.   SOLUTION  BY  THE  DISCRETE  MAXIMUM  PRINCIPLE  METHOD 

Let  the  demand  points  be  stages,  and  the  total  amount 
of  resource  which  has  been  transported  from  the  itn  depot 
to  the  first  n  stages  (demand  points)  be  state  variables 
x1?,  i  =  1,  2,  then 


l 


tn  „  xn-l  +  en    o  .  Q         4  =  w 
i    i      i'   l    '   l    i' 


(1) 


i  -  1,  2;   n  =  1,  2,  3,  4- 

It  must  be  noted  that,  though  there  are  3  depots  in  the 
problem,  there  are  only  2  state  variables.   This  is  because 
the  demand  of  each  stage  is  preassigned;  hence,  the  number 
of  the  units  supplied  from  the  3rd  depot  to  the  nth  stage 
can  be  obtained  by  subtracting  from  the  total  number  of 
units  required  by  the  nth  stage,  the  sum  of  the  units  sup- 
plied to  the  nth  stage  from  the  first  and  second  depot.   This 
is  equivalent  to  writing 

63  =  Dn  -  i    e?  ,  a-.l,  2,  3,  4.  (2) 

3  i=l 

Since  the  objective  is  to  minimize  the  total  cost  of  tran- 
sportation, we  define  this  objective  as  the  3rd  state  vari- 
able which  satisfies  the  following  performance  equation: 

xn=xn-l  +  I     Fn(en})   „  .  Q>  (3) 

J         J  i=l  J 

n  =  1,  2,  3,  4. 


It  can  be  shown  that  x_t  is  equal  to  the  total  cost  of  tran- 
sportation.  Hence.,  the  problem  of  minimising  the  total  cost 
of  transportation  becomes  that  of  minimizing  the  final  value 
of  the  3   state  variable,  xi,  by  the  proper  choice  of  the 
sequence  of '  8. ,  i  =  1,  2;  n  =  1,  2,  3,  4  for  the  process 
described  by  equations  (1)  and  (3). 

The  Hamiltonian  function  can  be  written  as 


Hn  =  I    z1?^-1  +  e")  +  z^u"-1  +  I    F"(en))  , 


n  »  1,  2,  3,  4.  (4) 

The  recursion  relation  for  the  components  of  the  adjoint 
vector  are  found  to  be 


,,n 


.n-1  _  _3H"  _.  _n 


ei     r£l~zi'    i  =  1'2  (5) 


1 
and 

ij'l,       n  =  1,  2,  3,  4  .  (6) 

Since  z.  and  x.~  are  considered  to  be  constants  in 
searching  for  a  stationary  point  or  in  finding  a  minimum 
point  of  the  Hamiltonian  function  given  by  equation  (4) 
with  respect  to  d^,    it  is  convenient  to  define  the  variable 
portion  of  the  Hamiltonian  function  as 


i=l        i=l  -1  •  1 


Since  the  performance  equations  (1)  and  (3)  are  linear  in  the 
state  variable  x?   and  the  coefficients  of  the  state  var- 
iables are  constants,  the  objective  function  is  absolutely 
maximum  (or  minimum)  if  and  only  if  H  is  absolutely  maximum 
(or  minimum).   Therefore,  the  optimal  sequence  of  the  de- 
cision vector,  6  ,  is  obtained  from  equation  (7)  if  and 
only  if  H  is  absolutely  maximum  (or  minimum) .   By  substi- 
tuting equation  (2-2)  and  (2)  into  equation  (7)  we  have 


nnnlQn 
2 


un   i   n  ,   n    n   „ni,ii|j    ,  n  ,   n    n   „,nnii 
Hv  =  \z-y   +  a^  -  a,  -  2b,D  )c),  +  (z„  +  a„  -  a,  -  2b, D  )Q 


+  (bj  +  b^)(e£)2  +  (bj  +  b^)(e£)2  +  a^Dn 


,  n  ,„n.  2    ,  n  n  n 

+  b  (D  )  -i-  2b3e1e2  . 

But  in  our  case  b  =  0  for  n  =  1,  2,  3,4;  therefore, 


,.n   ,  n  ,   n    n.Qn  ,  #  n  .   n    n, 
v  =  'zl   al  "  a3'  1   'z2    2  "  &V 


n,  „n 


,  n  .  nv2    n,  n ,2    n  n 
+  b1(61)   +  b2(82)   +  a3D   ,  n  =  1,  2,  3,  4  .    (&) 


Procedure  for  Optimal  Allocation : 


Step  1.   Calculate  the  breaking  value  zn  . 

i 


Stage  1: 

The  variable  part  of  the  Hamiltonian  equation  for "the 
first  demand  point  is 

h^  =  (zj-  +  0.5)0];  +  (z*  -  i)e^  -  .03(e^)2  +  100. 

Therefore  the  breaking  value  of  z^-  is  -0.5.   Taking  the 
partial  derivative  of  H:~  with  respect  to  8,  and  ea.uating 
it  to  zero  results  in 

3H1        , 

— -  =  0  =  zi  -   1  -  0.06  ei  . 

ae|  2  2 

However,  the  second  derivative  of  H  with  respect  to  9r 

v  2 

yields, 


& 


3(e2, 


-  -  0.06<  0 


3H1 
Therefore,  from  the  condition  — -  =   0,  we  cannot  obtain 

361 

1  2      1 

a  value  for  &2   which  yields  the  minimum  of  H  .   The  minimum 

of  H  occurs  at  the  boundary  of  the  constraint  of  0^,  as 
shown  in  Fig.  1.  Therefore,  the  conditions  for  H  to  be 
minimum  are 


H^;  -  min.  at:  Q*   =  0  if  2p2.j, 


02  =  50  if  z^<2.5. 


10 


31 


Thus,  the  so  called  breaking  values  of  z^  are  -0.5  and  2.5. 
The  six  conditions  at  which  Hy  may  be  minimum  are  presented 
in  Table  2. 

Stage  2: 

The  variable  part  of  the  Hamiltonian  equation  for  the 
second  demand  point  (stage)  is 

H2  =  (z{   -  2.5)6^  +  (z|  -  3.5)4  -  0.05(6?)2  +  350  . 

Therefore,  the  breaking  value  of  z|  is  3-5.   Taking  the 

2  ? 

partial  derivative  of  Hy  with  respect  to  9^   and  equating  it 

to  zero  results  in 

on  q  o 

— I  =  o  =  z    -  2.5  -  o.i  e^  . 

2  2 

However,  the  second  derivative  of  Hv  with  respect  to  6  yields 

a2Hv 

~   =  -0.K  0  . 

3(e2)2 

3Hv  2 

Therefore,  from  the  condition  — ^  =  0,  we  cannot  obtain  6-^ 

3d( 

o  2 

which  yields  the  minimum  of  H  .   The  minimum  of  H  occurs  at 
'  v  v 

2 
the  boundary  of  the  constraint  of  9-,  ,  with  curve  similar  to 

p 

Fig.  1.   Therefore,  the  conditions  of  Hy  to  be  minimum  are 


Table  2.   Conditions  necessary  for  Hv  to  be  minimu 


.1  :-• 


n 

Minimum  of  H  occurs  at 

V 

61 

82 

*1 

■§ 

0 

0 

>  -  0.5 

>  2.5 

o  s  e?-  <  50 

0 

=  -  0.5 

>  2.5  . 

1 

50 

0 

<  -  0.5 

>  2.5 

0 

50 

>  -  0.5 

<  2.5 

0 

50 

=  -  0.5 

<  2.5 

e|  =  0  or  50 

6^  =  0  or  50 

<  -  0.5 

<  2.5 

0 

0 

>  6 

>  3-5 

70 

0 

<  6 

>  3.5 

2 

0 

0  <  62  <  70 

>  6 

-  3-5 

70 

0 

<  6 

-  3.5 

0 

70 

1     6 

<  3.5 

6^  -  0  or  70 

e|  =  0  or  70 

<  6 

<  3.5 

0 

0 

>  -  2.5 

i-  2 

o  <  e^  <  ioo 

0 

=  -  2.5 

>-   Z 

100 

0 

<  -  2.5 

>  -   2 

0 

-  25  7?2   -    50 

>  -  2.5 

-6<232<_ 

2 

3 

o  <  e^  <  ioo 

-  25  z3  -  50 

=  -  2.5 

-  6  <7.\    <   - 

2 

o  <  q\   <  ioo 

-  25  z\   -    50 

<  -  2.5 

-  6<4   <- 

2 

0 

100 

>  -  2.5 

<  -  6 

0 

100 

=  -  2.5 

<  -  6 

9n  =  0  or  100 

1 

Q\   =  0  or  100 

<  -  2.5 

<  -  6 

13 


Table  2.   Conditions  necessary  for  Hv  to  be  minimum  (cont.) 


n 

Minimum 

of  Hy 

occurs  at 

e? 

eg 

n 
zl 

z2 

0 

0 

>  -  1 

>  1 

-  16.7  z^  -  16.7 

0 

-  5.8<ZJ<-! 

>  1 

80 

0 

<  -  5.S 

>  1 

0 

-0.5  z^  + 

i 

-  1  • 

-159£z^<2 

4 

-  16.7  z^  -  16.7 

-0.5  zj*  + 

1 

-5.8<z^<  -1 

-159  <z4  <2 
-159<z^<2 

0<6^  £  So 

-0.5  z^  + 

1 

<-  5.8 

0 

80 

>-  1 

<-  159 

-  16.7  z£  -  16.7 

O<02<  80 

-5.8<z^<-l 

<-  159 

ie£  =  0  or  80 

8^  =  0  or 

80 

<-  5.8 

<-  159 

14 


H2  =  min.      at:      92  =  0     if  z2 1  6    , 
e2  =  70  if  z2£6   . 

1  1 

Thus,  the  breaking  values  of  the  z?  are  6  and  3.5.   The  six 

o 
conditions  at  which  H  may  be  minimum  are  presented  in  Table 


Stage  3: 

The  variable  portion  part  of  the  Hamiltonian  equation 
for  the  third  demand  point  (stage)  is 

Hv  =  (zl  +  2-5'9i  +  <z2  +  2>92  +  -°2(e^)2  +  100  . 

3 
The  breaking  value  of  z:j  is  -2.5.   Taking  the  partial  de- 

3  3 

rivative  of  H^  with  respect  to  6£  and  equating  it  to  zero 

results  in 
3K.3 


— r  =  z\  +   2  +  .04  el 

de?   2        2 


e3  -  -25Z23  -  50  .  (9) 

Therefore,  the  upper  and  lower  breaking  values  of  z|  "are 

upper  t?2  =   -2,  when      0^  =  0, 

and  lower  z£  «=  -6,  when      9^  =  100  . 


15 


3  3       3  ? 

Hence,  H^  is  minimum  at  &2  ~  -25zi  -   50  for  -  6  zi     -2  as 

shown  in  Fig.  2.  The  nine  conditions  at  which  H-'  may  be 

minimum  are  presented  in  Table  2. 

Stage  4: 

The  variable  part  of  the  Hamiltonian  equation  for  the 
third  demand  point  (stage)  is 

h^  =  (z£  +  i)e£  +  (z1*  -  i)e£  +  o.o3(e^)2+  (e^)2  +  240 

Taking  the  partial  derivative  of  Hrj;  with  respect  to  9y  and 
equating  it  to  zero  results  in 

ngT  =  o  =  ztj-  +  l  +  o.o6  ej 


6^  =  -  16.7  z^  -  l6-7  •  (10) 


Therefore,  the  upper  and  lower  breaking  values  of  z  are 

upper  B7  ■=  -  1,    when     6^  =  0  , 

and  lower  z^  =  -  5.8,  when     0^  =  SO. 

Hence,  H^  is  minimum  at  87  -  -  16. 7  z^  -  16. 7  for  -  5.8  < 
ztj-<-  1  with  curve  similar  to  Fig.  2.   Taking  the  partial 
derivative  of  H^:  with  respect  to  9^  and  equating  it  to  zero 
results  in 
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Fig.  2.   H^  v.s.  e\. 
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9H4      ■         i 
30^       2        2 


92  =  "  °'5  Z2  +  °-5  -  (11) 


Therefore,  the  upper  and  lower  breaking  values  of  zi,   are 

upper  zT,  ~   1,     when     Gg  =  0, 
and  lower  z^  =  -159,   when     9^  =  80. 

Hence  H^  is  minimum  at  6^  -  -  0.5  z'j;  for  -159 -Zg-1  with 

curve  similar  to  Fig.  2.   The  nine  conditions  at  which  H 

v 

may  be  minimum  is  given  in  Table  2.   The  conditions  for  all 
Hy  to  be  minimized  are  summarized  in  Table  2.   The  breaking 
values  of  all  the  stages  are  summarized  in  Fig.  3a  and 
Fig.  3b. 


Step  2.   Systematically  searching 

Each  combination  of  the  interior  and/or  the  breaking 
values  of  z-,  and  z,  are  systematically  searched  for  a 

feasible  solution;  cases  which  do  not  satisfy  the  con- 

*■   n 
straint  £     6-,  =  Wi  will  be  eliminated. 

n=l 


a)   Region  to  be  searched 

In  this  searching  procedure  regions  or  points  of  the 
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Fig.  3a.   The  breaking  value  of  adjoint  variable  z-. 


n=4 


4     12 


1   2.5  3.5 


Fig.  3b.   The  breaking  value  of  adjoint  variable  z~. 
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breaking  value  of  z?  are  searched  systematically  by  exhaust- 
ing all  possible  combinations  in  the  region  of  the  breaking 
value . 

b)   Elimination  of  infeasible  solution 

Most  of  the  infeasible  solutions  result  from  not  sat- 
isfying the  constraint  we  put  into  the  system,  i.e., 

'<   n 

I    e1  =  w   i  =  1,2. 

n=l  x    1 

Example  1.   The  combination  of  regions  z^  and  of  zg  given  by 
-  0.5  <z1<   6;  2.5  <k2<  3-5: 

From  Table  2,  the  tables  of  conditions  necessary  for 

H"  to  be  minimum,  we  have  the  conditions  shown  in  Table  3. 

£   n 

It  is  obvious  that  this  is  an  infeasible  case,  since  £  9,  = 

n=l  J 

50  +  0  +  100  +  80  =  230  is  larger  than  W, ,  where  W  =  80. 

However,  there  are  cases  where  the  value  z-  calculated 

from  the  equation  of  optimality  (such  as  equations  (9),  (10), 

and  (11)  when  the  value  of  8*    are  known)  does  not  fall  in 

the  region  of  search. 


Examo 


le  2.   The  combination  of  region  for  z  and  z2  given  by 


-  2.5  <  z  <  -  1;  -  6 <  z2  <  -  2: 

From  Table  2  we  obtain  the  conditions  necessary  for  Hv 

to  be  minimum  as  shown  in  Table  4-   The  condition  of  6-j_  =  0  or 

50  and  8.  =  0  or  50  means  that  the  minimum  value  of  H  is 
2  v 

either  at  e|  =  0  and  el,   -  50  or  8^  =  0  and  6*  -  50. 


20 


Table  3.  9i  corresponding  to  the  condition  that  z,  and 
z"  are  in  the  region  of  -0.5  <z    16  and 


2.5 


<«n 


<3-5 


\    i 
n  X. 

1 

2 

3 

Dn 

1 
2 
3 
A 

0 
6^  =  0  or  70 
0 
0 

0 

e|  »  0  or  70 

0 

0 

50 

0 

100 

80 

50 

70 

100 

80 

Wi 

130 

90 

80 

300 

Table  k  •      6.  corresponding  to  -2.5<z?<-l  and  -6lz2<-2 


Xs^   i 
n  \^ 

1 

2 

3 

Dn 

1 
2 
3 
4 

s|  =  0  or  50 
9f  -  0  or  70 

0 
-16. 7z^  -  16.7 

62  -  0  or  50 
BJ   «"  0  or  70 
-25  z^  -  50 
-0.5z|  +  1 

50 

70 

100 

80 

Wi 

130 

90 

80 

300 
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By  using  equation  (3)  we  can  compare  the  value  of 


Hv  for  both  cases  as  follows: 
For  G?-  =  0,  0^  =  50 


1  „1 


upper   Hy  =  {z\  +  &\  -  z\)q\  +   b^Og)2  +  &\   D 

=  (-2.1  +  1  -  2)  50  -  .03  (50)2  +  2  x  50 

-  3.1  x  50  -  75  +  100 
=  -  130 

lower   H1  =  (-5.9  +  1  -  2)50  -  .03(50)2  +  2  x  50 
=  -  6.0  x  50  -  0.03  x  (50)2  +  2  x  50 

-  -  345  -  75  +  100 

-  -  370 

For  9*  =  50,  q\   =  ° 


upper  Uy   -  (-  1.1  +  2.5  -  2)  50  +  2  x  50 
-  -0.6  x  50  +  100 
«  -  30  +  100 
=  +  70 

lower   H*  -  (-  2.4  +  2.5  -  2)  50  +  2  x  50 
=  -  1.9  x  50  +  2  x  50 
=  +  10 

It  is  obvious  that  6^   =  0,  92  =  50  gives  a  minimum,  since  the 
value  Hv  is  negative  throughout  the  regions  of  z2  while  the 
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value  of  Hv  when  8..  =50,  82  =  0  is  positive  throughout  the 

region  of  z,.   Similarly,  we  assign  8^  =  70,  02  »  0;  therefore, 

9!*   =  W,  -   J"  p"=  130  -  70  =  60.   By  using  the  equation  of 

1    1    n=l  X 
optimality  given  by  equation  (10)  we  find  that 

-16.7  z£  -  16.7  =  60, 


z4  =  60  +  16. 7  _  .  4>6. 
1     -16.7 


By  equation  (6)  we  have 

z,  =  z-,  =  zf  =  -  4-6   for  n  =  1,  2,  3,4- 

Since  the  region  of  z-,  we  are  now  searching  is  -  2.5<z^<-l 
and  since  the  value  of  z.  we  just  calculated  is  out  of  the 
region,  this  is  an  infeasible  combination  of  the  region  of 
search. 

c)   Feasible  solution 

Feasible  solutions  are  combinations  of  regions  of 
z^  which  satisfy  all  the  constraints  of  the  system  and  the 
value  of  z?  calculated  from  the  equation  of  optimality  is 
within  the  region. 

Example  3 •   A  feasible  solution  corresponding  to  the  values 
of  z"  at  the  point  of  z"  «=  -  2.5  and  z£  in  the  region  of 
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6<zn<-  2  which  satisfies  the  constraints  is  presented 


in  Table  5 
By  co: 


By  comparing  H  as  we  did  in  example  2  we  have  0-^0, 


62  =  50;  similarly,  we  have  6;L  =  70,  92 


0. 


Since  z"  -  z^  ■=  -  2.5,  from  equation  (10)  we  have 

64  =  -  16.7  z*l   -  16.7 

=  -  16.7  x  (-  2.5)  -  16.7 
-  41.7  -  16.7  =  25 

ei  =  wi  -  el  -  ei  -  ei 

=  130-0-70-25 
-  35. 

From  the  constraint  W2  -  90  and  z2  =  z|  ■  z4  we  obtain 
50  +  0  -  25  z2  -  50  -  0.5  z2  +  0.5  =  90 


or 


z=.  89^5  =.3. 51 
*  25-5 

which  satisfies  the  condition,  -  6  5z2<-  2.  The  corres- 

ponding  value  of  Gs  obtained  by  substitution  of  z     =   -  3<51 

is 

Q\   =  "  25  (-3-51)  -  50 
=  87.6  -  50  =  37.6 
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Table  5.   6^  corresponding  to  the  condition  that  z"  and  z" 
are  in  the  region  of  z"  ■=  -2.5  and  -6 £z™  <-2. 


n   \ 

1 

2 

3 

Dn 

1 

e|  =  0  or  50 

62  =  0  or  50 

50 

2 

e|  =  0  or  70 

e|  =  0  or  70 

70 

3 

o  le^lioo 

-  2.5  z^  -  50 

100 

4 

-16.7  z£  -  16.7 

-  0.5  z1*-  +   0.5 
2 

80 

»i 

130 

90 

go 

300 

Table  6.   6^  corresponding  to  z"  =  -  2.5  and  -6  i  z?  < -2 


1 

2 

3 

Dn 

1 

0 

50 

0 

50 

2 

70 

0 

0 

70 

3 

35 

38 

27 

100 

4 

25 

2 

53 

80 

W. 

130 

90 

80 

300 

25 


or 


el  -  38 

82  =  "  °*5  (-3-51)  +  °-5  =  2.25 


or 


-2-  2- 

The  value  of  9.,  can  be  found  by  using  the  constraint, 

t    e"  =  d". 

i=l 

The  itemized  cost  for  the  above  solution  is  presented 
in  Table  7,  and  the  total  cost  is 

4   3   n  n 

E   E  F?(8?)  ^  171.25  +  125.8  +  186  -  $483.05. 
n=l  1=1    x 

Another  feasible  solution  corresponding  to  the  values  of 
z"  in  the  region  of  -  5-8<zn<  -  2.5  and  z5  in  the  region 
-  6  5zi<-  2  which  satisfies  the  constraints  is  presented 
in  Table  8. 

The  allocation  of  the  value  of  6i  and  the  calculation 
of  transportation  cost,  similar  to  example  3,  are  shown  in 
Tables  9  and  10. 

The  total  transporation  cost  for  this  combination  of 
regions  is 


3    4 

,£,  t  Me")  =  $173.30  +  125.80  +  184.00  -  $483.10. 
i-l  n=l 
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Table  7.  Cost  of  transportation  corresponding  to  z^  ■  -2.5 
-  6iz£<-  2 


\^                     i 
n          >v 

1 

2                                  3 

nQn 
al°l 

,  n,„n.2 

b1(e1) 

nQn 
a262 

,  n,_n,2 

n„n 
a->6n 

n,   n.2 

b3(e3) 

1 
2 

3 
4 

0 
175 
122.5 
100 

0 
-245 
0 
18.75 

50 

0 

114 

4 

-75 
0 
28.8 

4 

0 

0 

27 

159 

0 
0 
0 
0 

397.5 

-226.25 

168 

-42.2 

186 

0 

n^l 

171.25 

125.8 

186 



n=l  i=l         ■"■ 

483-05 

Table  8.  9^  corresponding  to  the  condition  that  z,  is  in  the 
region  of  -5.8  5.z-,  <-2.5  and  z,  is  in  the  region  of 


-  6<Z2  <-2 

>v     i 

n  \. 

1 

2 

3 

Dn 

1 

ej  =  0  or  50 

62  =  0   or    50 

50 

2 

ef  =  0  or  70 

e|  =  0  or  70 

70 

3 

o<e^5ioo 

-    25   z|   -    50 

100 

4 

-  16.7  z£  -16.7 

-  0.5  zjj  +  0.5 

80 

Wi 

130 

90 

80 

300 

27 


Tabl 


e  9.  6-  corresponding  to  -5.8  <a"  <  -2.5  and  -6<z2<-2 


Nu  i 

1 

2 

3 

Dn 

n   N. 

1 

0 

50 

0 

50 

2 

70 

0 

0 

70 

3 

34 

33 

28 

100 

4 

26 

2 

52 

80 

W. 
i 

130 

90 

80 

300 

Table  10.  Cost  of  transportation  corresponding  to  -5.8  <z 
<  -2.5  and  -6iz£<  -2 


n^\^ 

1              2              3 

nQn 

a1e1 

.  n/Qn \ 2 
bj  (e1) 

nQn 
a262 

,  n,  n,2 
b2(82) 

afef 

, n,„n>2 
b^O,) 

1 
2 
3 
4 

0 
175 
119 
104 

0 
-245 
0 
20.3 

50 

0 

114 

4 

-75 
0 
28.8 
4 

0 

0 
28 

156 

0 
0 
0 
0 

398 

-224.7 

168 

-42.2 

184 

I  F^e-1) 

n=l 

173.3 

i 

125.8 

I84 

I  £#* 

483.10 

2fi 


Step  3-  The  optimal  solution  is  obtained  by  comparing  the 

total  transportation  cost  of  all  the  feasible  regions.  In 

our  case  there  are  only  two  feasible  regions.   Since  the 

total  transportation  cost  corresponding  to  z?  =  -2.5  and 

-6iz"<-2  is  less  than  that  of  -5.gz?<-2.5  and  -6<z?<-2, 
2  **  * 

the  optimal  solution  is  &x   corresponding  to  z-^  =   -2.5  and 
-6iz5<-2  which  is  shown  in  Table  6,  and  the  minimum  total 
transportation  cost  -  $4&3-05. 
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h.      COMPUTING  PROCEDURE 

Step  1.  Read  in  data:   the  values  of  ai,  bn,  W^,  and  D  . 
Step  2.  Calculate  the  breaking  value  of  zn. 

The  equation  used  for  the  calculation  depends  on  the 
value  of  b?  as  follows : 

n  n 

a)  For  b^<0,  the  secondary  derivative  of  H  with 

respect  to  0^  will  be  negative.   Therefore,  from  the  con- 
dition  — ~  =  0,  we  cannot  obtain  a  9.  which  yields  the  mini- 

n  n 

mum  of  H  .   The  minimum  of  Hv  occurs  at  the  boundary  of  the 

n 
constraint  of  Q±   as  shown  in  Fig.  1;  the  equation  to  be  used 

in  this  case  is  as  follows: 

n 
zi  ~  ~  ai  "*  as  "  2  bi  ~2  (12) 


b)  For  b.  =  0,  the  cost  function  is  linear,  and  the 
minimum  will  occur  at  the  boundary  of  the  constraint  of  8^. 
The  breaking  value  is  given  by 

4   -  -  4   +  a£  .  (13) 

c)  For  bj.  >  0,  the  cost  of  transportation  is  non-linear 
and  the  second  derivative  of  H"  with  respect  to  9i  is  larger 
than  zero.  The  shape  of  the  curve  of  Hy  verses  9?  is  as  shown 
in  Fig.  2.   The  minimum  of  Hy  can  be  obtained  by  equating  the 
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derivative  of  h"  with  respect  to  e±   to  zero.      The  lower 
breaking  value  (Lower  z")  and  the  upper  breaking  value 
(Upper  z?)  are  at  the  boundary  of  the  constraint,  that  is, 
6?  =  Dn  and  6?  =  0  respectively,  or 

Lower  zn  =  -  a1?  +  an  -  2bn  Dn  ,  (H) 

i      1    s     i     >  K±Ll1 

Upper  z1}   =  -  an  +  an  .  (15) 

Step  3.   Define  the  regions  of  search  according  to  the 
breaking  values. 

This  step  consists  of  three  substeps. 

a)  Rearrange  the  upper  breaking  value,  z?,  in  an  order 
of  largeness  and  represent  it  by  ZEU^',  where  J   =  1,  2,  ..., 
indicate  the  order  of  largeness.   The  lower  breaking  value 
of  the  nfc   stage  will  use  the  same  j  value  as  the  upper 
breaking  value  of  that  stage  and  represent  it  by  ZBlA 

b)  Assign  the  upper  limiting  value  of  regions  to  be 

searched.   The  upper  limits  of  regions  to  be  searched  are 

assigned  according  to  the  value  of  bb?  (the  bn  value  in  the 

l       l 

new  index)  as  follows: 

(1)  bb^<0,  ZOP^  =  ZBU^'  -  0.005  (16) 

(2)  bb^'  =  0,  ZOP^  »  ZBU?  .  (17) 

Z0P^+1  -  ZBU^  -  0.005  (18) 
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(3)      bb?>0     Z0PK  -  ZBU^  -  0.005  (19) 


Z0P^+1  =  ZBL^   -  0.005  (20) 

The  value  0.005  is  a  small  increment  between  the  lower 
limit  of  the  previous  region  and  the  upper  limit  of  the  cur- 
rent region. 

In  the  search  for  a  feasible  solution,  the  upper  limit 
of  each  region  will  serve  temporarily  as  the  zj  value  in 
the  first  trial. 

k,  the  new  index  for  the  regions,  follows  the  order  of 
the  j  index.   Therefore,  ZOPJf,  k  -   1,  2,  ...,  7,  is  not  in 
order  according  to  the  largeness  of  the  value  of  ZOPf. 

c)  Rearrange  ZOP^,  k  ■  1,  2,  ...,  7,  in  an  order  of 
largeness  of  the  value  of  ZOPl?. 

The  results  of  step  3  are  shown  in  Figs.  4a  and  4b. 

Step  4.   The  first  allocation  of  e". 

The  optimal  value  of  9n  which  corresponds  to  the  value 
of  ?,±,   given  (ZOP^),  i  =  1,  2,  can  be  classified  into  the 
following  5  cases  (see  Table  2) 

„k 


a)  e"  =  0,     ZOP*  >z?  . 

b)  6?  =  -  (ZOP*  +  a£  -  a£)/(2  bj)  , 

Lower  z?l  Z0P^<  Upper   z"   . 


b>0 
3-3 

n=4 


k   5- 


j.jo  7- 


-5.8 


b=>0  b>0  b=0 
4     3  2 


3     4  1 
r6  I  3 1 

4  J 

'-5  I 


' 


t 


I  I 


-2.5  -1  -0.5 
n 


Fig.  4a.   The  region  or  point  of  search  of  z. . 
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b<0 

1  Step  3a 

2  . 
! 


Step  3b 
Step  3c 


k=5  " 


!,&•?' 


b>0   b<0  b=0 

3  2   1  Step  3a 

4  12 


-  3-i 


-Z-j 


-  3 


-2-! 


Step  3b 
Step  3c 


2.5  3.5 


n 
Z2 

Fig.  4b.   The  region  or  point  of  search  of  z2. 
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c)  e"  =  Dn,    zopIc<  zn 


d)  O£0i<D  ,   ZOP.  =  zi.   In  computer  programming 
the  optimal  value  6^  is  represented  by  Dn  +  1,  that  is, 

6.  =  Dn  +  1  at  this  step. 

e)  e?  =  Dn-[-(Z0P^a  +  4a-a£)/2b?a]    , 


ZOP/x  z"  and 
i      i 

Lower  zn  <  ZOP.  <  Upper  z.   .   In  computer  programming  the 
optimal  value  0n  is  represented  by  Dn  -f  2,  that  is  e"  =  Dn  +  2 
at  this  step. 

Step  5.   Preliminary  elimination  of  infeasible  combinations 
of  regions  of  search. 

Most  of  the  combinations  of  the  region  of  search  are 
eliniminated  in  this  step.  This  is  done  by  using  the  op- 
timal values  of  8.  which  have  been  allocated  in  Step  4.1a 
and  4.1c  (8j  =  0,  e"  =  Dn),  and  the  constraint 

N 
(W^  =  £  6^ ) .  The  two  conditions  that  give  the  infeasible 
n=l 

regions  are 

N  n    n 

a)   £  e^>ViT^;  this  may  occur  when  we  have  6j  =  D 
n=l 

for  several  n's  in  the  ith  origin,  and 
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i   $   n  n 

b)   2.  ei<wi>  this  may  occur  when  v;e  have  8^  =  0  for 
n=l 

many  n's  in  the  i   origin,  even  though  we  assign  the  maxi- 
mum value  (Dn)  for  the  rest  of  n's. 

Step  6.   The  second  allocation  of  e":   Dn  +  1,  Dn  +  2, 

6n  =  -  (ZOPJ  -i  a?  -  an)/2bn  . 
-1  1     1     s     1 

This  step  is  divided  into  7  substeps. 

1)  The  identification  of  the  different  conditions  of 
optima lity. 

In  each  origin  the  number  of  occurences  of  each  type  of 
optimal  condition  classified  in  Step  4  is  counted.   The  type 
of  optimal  condition  at  each  stage  is  also  identified.  The 
counters  and  indexes  for  different  types  of  conditions  of 
optimality  are  shown  in  Table  11. 

2)  Find  the  available  resource  for  the  unassigned  stages. 
The  difference  between  the  total  resources  available 

at  each  origin  or  depot  (W^)  and  the  sum  of  all  0?  =  Dn 
assigned  (SUM  D^ )  is  the  available  resource  for  the  unas- 
signed stages  (9  SUM) .   It  can  be  put  into  computer  equations 
as  follows : 

SUM  D±   =  SUM  DA  +  Dn  (21) 

9  SUM  =  Wj  -  SUM  D^^  (22)  ' 


Table  11.   Indicator  in  computer  program 
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Type  of 

Con- 
dition 

Kind  of 
Condition 

Temporary 
Value 

No.  of 
Types  at 
each  Origin 

Sum  of 
2 

J   ke  etc. 
i=l 

Index  to 
identify 
stage 

4.1b 

Equation  of 
optima lity 

0<  en<  Dn 

ke 

ckk 

meke 

4.1e 

Severely 
Constrained 

e?  «  Dn+2 

kc.kk 

P 

mkkc 

4.  Id 

Free  varied 
of  Theta 

e?  =  Dn+i 

k 

q 

mk 
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When  the  assigned  stages  include  the  temporary  value 
of  8n  »  equation  of  optimality  and  9^  =  Dn,  the  sum  of  the 
resource  of  the  assigned  stages  (SUM^)  will  be  the  sum  of  all 
6n  -  Dn  and  9n  =  -  (z?  +  a?  ~  a")/(2b^).   The  available  re- 
source  for  the  unassigned  stages  (TH^)  will  be  the  difference 
between  the  total  resource  available  at  that  origin  (W^)  and 
the  sum  of  the  resource  of  the  assigned  states  (SUK^),  which 
can  be  written  in  computer  equations  as  follows: 

sm±   =  SUI^  +  Dn  (23) 

SUJ^  -  SDMjl  +  b"  (24) 

THi  =  W±  -  Smi±  (25) 

When  the  assigned  stages  include  the  temporary  8*  value 
(09-^)  of  the  severely  constrained  condition  and  0^  =  Dn,  the 
available  resource  for  the  unassigned  stages  (6S^)  will  be 
the  difference  between  the  total  resource  available  at  that 
origin  (W. )  and  the  sum  of  the  resource  of  the  assigned 
stages  (SUM  06.).   It  can  be  written  in  computor  equations 
as  follows: 

sum  6d±  =  sum  ee±  +  Dn  (26) 

SUM'SSj  =  SUM  &dx   +  90n  (27) 

esi  =  wA  -  sum  edi  (28) 
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3)   Allocation  for  6^  which  is  given  by  the  equation  of 
optimality 

6?  =  -  (*?  +  a"  -  an)/2bn  .  (29) 

1  l     1-    8     i 

The  z^  value  used  in  finding  the  optimal  G.   may  be 
divided  into  5  catagories. 

a)  At  the  point  of  search.  When  there  is  a  condition 

of  free  variation  of  6-  in  the  origin,  the  ZUP.  should  be 
i  i 

used  immediately  and  is  fixed. 

b)  In  the  region  of  search.   In  this  case  we  have  to 
search  a  whole  region  or  a  particular  interval  of  zn  value 
to  find  the  optimal  value  of  z . .  The  equations  for  finding 
the  optimal  z.  depend  upon  the  number  of  equations  of  opti- 
mality in  the  same  origin  as  follows: 

(1)  One  equation  of  optimality  condition  (ke  =1). 
.im  ,  /„im   ,  im 


z. 


+  la-   -  a,  )  N 


en  =,  _i f  "  3   =  _  (w  _  j-  Dnj  K   _  QSm     ^Q) 

1       2  b\m  n=l 

n^im 


2  b*m  6SUM  -  (anim  -  a™1) 

1  AS 


(2)  Two  equations  of  optimality  conditions  (ke  =  2). 

Z.in+  (a*m  -  a*m)   zf1  +  (a>  -  a^)  N 

-i i— 2 h  _i 3 § =  _  (W.    _   V   Dn) 

2  bjm  2  bf  (  X       n=l 

im'im 
n^jm 

=■  -  esuM  , 
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irr.  im      ,  im.  im  im,  ,im,  jm  im, 
2b,  b?  eSUM  -b,  (a  -a=  )-b,  (ar  -a^  ) 


z  „zim=  zjm  „  -"i  "l     .11   s    l   l   s 
1    i     *  (b^W1) 


(3D 
(3)   Three  equations  of  optimality  conditions  (ke  =3)' 


«*m  +  (a*m  -  aim)   z J™  +  (aJm-aJ'm)   zkm  +  (a^-a1™ 
l      i s   +  l      l   s   +  _i _i s_ 

2  bf  2  bf  2  b*m 


(w.  -    y   d11       )  =  -esuM    , 

1   n-l 
•  -    nf-im,  jm,km 


or 

„,  inujin.  km_om,  ,im,km,  im  im, 
im  _  jm   _  km  _  -2bi  Pi  *>i  6SUH  -bj  bj  (&1  -as  ) 
zi   zi    zi   '"  zi  bimbJm+bombkm 


,im,km,  im  jm,  ,im,  jm,  km  km. 
-bj  b±    {a(   -aJs    )-b  bJ  (3i  -as  ) 


..  inukra 

+bi  bi 


£JL  .   (32) 


(4)  Four  equations  of  optimality  conditions  (ke  =  4)'. 
zf+faf-a^)   Zf(af-af)   ■{B+(aJm-«JmJ 


+ 


2  b|m         2  b^'m  2  b^1 
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+  "      -   ■  1m  =  "    (Wi   -     I     D    >   "    ~eSUW     ' 


2  bf  n=l 


n^im, jm,km,lm 


3fm  -  zjm  =  zkm  =     lm 
iii2 


=  -(2bimb^bim6SUK     +  bJBbJmbJBl(aJ»-aJm) 

+  bfbfb^af-af )   +  bfbf  b^faf-a^) 

+  bjmbfb^(aKalm))/     b^bfb^m  +  bimomblm 

X         X         X  X  lb  XXX  XXX 

+  birabknlblra  +  bJ'mbkrablm   .  (33) 

ill     iii  \jji 

Similarly  any  number  of  equations  of  optimality  condi- 
tion in  the  same  origin  can  be  derived. 

c)  Upper  limiting  value  in  the  region  of  search.   When 
there  is  a  severely  constrained  condition  (4.e)  in  the 
origin  and  SS^O,  the  maximum  value  of  z,  in  the  interval 

of  the  region  that  we  are  searching  is  used.   This  value 
corresponds  to  the  minimum  value  of  9.  that  can  be  allocated 
to  stage  n  in  this  region  of  search.   This  occurs  when  9S1 
from  equation  (28)  is  less  than  or  equal  to  zero. 

d)  The  region  of  search  with  severe  constraint.   When 
there  is  a  severely  constrained  condition  in  the  origin  and 
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the  value  es^^  from  equations  (26)  is  larger  than  zero,  z± 
of  the  equation  of  optimality  (29)  can  be  calculated  from 
(6.3b)  instead  of  (6.3c)  and  the  value  zn  can  be  anywhere 
within  the  region. 

e)   Iterated  value  of  z".   When  there  is  more  than  one 
origin  with  the  property  of  (6.3d)  we  have  a  different  series 
of  values  of  9?  each  time  we  repeat  the  DO-loop  of  step  6. 
This  is  because  the  value  z.  is  changing  each  time  we  com- 
plete the  DO-loop  of  step  6.   Therefore  the  DO-loop  of 
step  6  has  to  be  repeated  until  we  get  a  set  of  constants, 
z",  which  give  us  the  optimal  6i  value  for  this  combination 
of  regions  of  search.   When  this  is  reached  the  triggering 
counter  (ckk)  will  decrease  to  zero  and  the  next  substep 
will  be  executed. 

4)  Allocation  for  6.  with  severely  constrained  con- 
dition case. 

The  reason  this  case  has  a  severely  constrained  condi- 
tion is  that  the  condition  for  optimal  Q^   is  0i  =  D  ,  but 
it  is  constrained  by  the  priority  we  give  to  the  equation 
of  optimality  in  the  preceding  or  the  next  origin.   If  this 
condition  is  nowhere  present  in  the  total  system  at  the 
current  combination  of  regions  of  search  (p=0),  this  sub- 
step  will  automatically  skip  to  the  next  substep. 

The  maximum  allowable,  available  resource  is  allocated 
to  6j  in  this  case  so  that  the  value  of  9^  is  as  close  to 


w 


the  opti.mal  condition  as  possible.  This  is  done  by  checking 
both  constraints 

imk  „  Dn  _  n  (35) 

1  13 


and 


efrak  =  w.  -  ?  en  .   (36) 

i     i   A   i 

n-^irak 


and  choosing  the  smaller  of  these  two  values.  This  value 
is  the  maximum  allowable,  available  amount  of  resource. 

5)   Allocation  of  8?  with  free  variation  of  6.  condition. 

This  is  the  last  type  of  9^  in  the  active  origin  to  be 
allocated.   Since  the  optimal  value  can  be  anything  from 
0  to  D  ,  it  also  serves  as  an  allocation  to  satisfy  the 
constraint 

N    n 


When  this  property  is  not  present  at  the  region  of  search 
(q=0),  this  substep  will  automatically  skip  to  the  next 
substep. 

6)  Control  of  Allocation.  The  control  of  the  alloca- 
tion serves  two  functions: 
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(1)  It  serves  to  check  if  the  optimal  z±   is  within  the 
interval  of  the  region  of  search,  and 

(2)  as  a  trigger  mechanism  for  substeps  (6.3),  (6.4), 
and  (6.5)  which  can  be  described  as  follows: 

For  substep  (6.3)  when  the  counter  is  triggered, 
(ckk=0)  the  search  for  optimal  O7}   is  complete 
for  this  substep  and  it  will  go  back  to  step  6 
again.   This  time  the  substep  (6.4)  is  executed  if 
p>0  but  if  p=0  the  next  substep  (6.5)  will  be 
executed  if  q > 0  and  if  q=0,  the  computer  will  go 
to  substep  (6.7) 
7)  Allocation  of  constrained  origin  (s  -depot) 
The  reason  we  called  this  a  constrained  origin  is  that 
the  number  of  the  units  supplied  from  the  constrained  origin 
(sth-depot)  to  the  nth  stage  can  be  obtained  by  substracting 
from  the  total  number  of  units  required  by  the  n**1  stage  the 
sum  of  the  units  supplied  to  the  nth  stage  from  the  first 
through  (s-1)   depot  (or  active  state).   This  is  equivalent 
to  writing 

e"  =  Dn  -  siX  en  =  d"  -  b?  -  e5  . 

s        ^  1         1    2 

Step  7.   Calculate  total  cost  of  transportation 
i 
The  cost  of  all  feasible  solutions  is  calculated  by 

substituting  the  optimal  9i  into  the  cost  function  equation, 
and  summing  all  the  costs  of  transportation  of  each  feasible 
solution. 


A3 


Step  8.  ■  Compare  the  total  cost  with  that  of  the  previous 
minimum. 

The  total  cost  of  transportation  of  the  new  feasible 
solution  is  compared  with  that  of  the  previous  one.  The 
new  minimum  is  then  stored  for  future  comparison. 

Step  9-  Punch  out  output  on  cards 

All  the  values  of  the  feasible  solution  are  punched  on 

n 
cards,  including  9^  i  =  1,  2,  3  and  n  =  1,  2,  3,  4;  the 

region  of  search  of  the  new  minimum;  the  up-dated  minimum 

cost;  and  the  current  total  cost. 
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APPENDIX  A  COMPUTER  FLOW  CHART 


The  steps  of  the  computer  program  are: 

1.  Read  in  a J,  b£,  aj»,  Dn,  w . ,  »3 . 

2.  Calculate  the  breaking  value,  2". 

3.  Order  the  breaking  values  from  the  largest  to  the 
smallest;  and  assign  the  region  of  search. 

4.  Choose  a  pair  of  the  region  of  search,  and  make  the  first 
allocation  of  e^1;  0,  Dn  +  1,  o"  -  -(ZOP?  +  a£  -  an)  ,  Dn 


2bi 


and  Dn  +  2. 


5.  Check  the  feasibility  of  the  region. 

6.  Make  the  second  allocation  for  8n:   en  =  Dn  4  1; 

Qn   -(ZOPJ  +  an  -  an)        n    n 

6H  =  i    i    s  ;   and  9  =  Dn  -l  2. 

n  x 

2b. 

7.  Calculate  total  cost  of  transportation. 

8.  Compare  the  total  cost  with  that  of  the  previous  minimum 
cost. 

9.  Punch  out  the  6i  which  gives  the  current  minimum,  where 
i  =  1,  2,  3  and  n  =  1,  2,  3,  4;  the  region  of  the  new 
minimum;  the  registered  minimum;  and  the  current  minimum. 
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Stpp  I   r.  2 


INITIALIZATION 
SHIM   =  1,000,000 

s  =  3 

IT  =  4 


C  DO   i  =  1^2) 


-J 


DO  n=  1,4 
IZZZUZZ 


> 


R2AD  a^b^a^D11,;^/^ 


<0 


«J  =  Dri/2 

UPPUR  Zn  =  -aV+bW 
1            1     s     :u  1 

liOVV.Ul  Z?  =  UPP3R  Z? 
1                       1 

UPPER  Z*  =  -a?  +  a* 


UPP3H  Z?  =  -a1.1  +  a11 


LOV/SR  Z"  =  UPP3H  ZV 


PUNCH  1.H.UPP3H  Z1?, LOWER  Z? 


COUTH:  U 


GO   50  ST2P  3 


k& 


Step  3 


(pO  3   =  1,4 


OT^  =  -1,000 


(j")0  n  =~T^}- 


SiU^-Upper  Z? 


>    0 


couTiirj]/ 


Ztf?°  =  zopno 


LTJp1?0  =  -1,000 


mLl     =  lower  £?° 


PUITGH  i, j.ZBOj 


CONTINUE 
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INITIALIZATION 
k  =  0 


"fk  =  k+1 
ZOP1?     =  ZB&| 


k+1 


ZOP^     =  Z31j|-.005 


k  =  k+1 
ZOPk     =  ZBU-!  -.005 


k     =  k+1 


20Pk     =  ZBI^-.005 


nk.=  k 


-   COHCTTOB 


t 


puijcii     nk. 


T 


C0KT1MDE 
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— (jX)      jo  =   1,    njo) 


HfflCH   ZOTJ1 


coittinue 


CO   TO  STEP   4 
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STJPJ 

INITIALIZATION 
0=0,    p  =   0,    q   =   0 
COUNT  =  o,    okk     =  0 


=  0,      ZZp  a  0 


SUE,    =  0,    SUK2   =   0 
SW.I3   =  0,    SUBIffj    =  1VT 

qTTTIW        -    "'  ^TTT^'f    —    "r 

OUlWp     -     .'rp,      .JU!       -7-     "  rji 

kkk  *■  0  PMJ  =  0 

i  =  1,2,3  n-1, 2,3,4 


PUNCH  jo,  j,  ZOPJ   ,  ZOP- 

x-i     i 


UPPER  Z^  =  OTj, 
UPPER  ^'•1   =  2U? 


<0 


i:;opv  -  upper 


>0 


<0 


<  0 


<C 


<?!?  =  Nn  ■- 1 

e£_1-  s  +  i 


q;  =  d"  +  l 


_^o   ( 52 


?JV^BR  I,. 


e£  =  d"  +  i 


aSFjfj-UPPEH 


M 


e£  =  o 


i-1 


=  D  +1 


0 


9?   =   0 

1 


<o/ i0  -7T\>0 


«g  ,,-(^!i^i-4? 


2tf 
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r 


<  °  /zr.      ~~z\  > c 

ZOPJ-LOWER   Z?) 


<r5-UcpJ°ruppER  zj_^c 


<  0 


zcp"J°-Lowsa  zn  ,)— i 
i-i  1-1/ 


\    >0 


e?  =  Dn  +  2 


2b?_l 


@     ' 


;r 


„     -(15pJ+a?-an) 

a"  1 i si 

1  2b$> 

8?-l  =  0 


_;        r~\  >o 

^Q-LOKSR   Zi.y J 


6i    = 


■(ZOPi+a^-a^) 
2b?  ~ 


e"  =  d"  +  2 


'i-1 


, J      0      n. 

-(ZOPfo^-Sg) 

2bn 

.(^OPJo  +an      -a"1 
2b"   , 


5Jc   .t„n      _«ninIU   Hn      mn%2 


Hl-1   -    (ZOP^+a^-a*)^   bJ_x(D 
Hi        =    (ZOPJ   +   4  -   4)Dn  +   b?(D")2 


H2j 


53 


en  =  o 
i 


PM 


1-1 


=   1 


% 


No 


8?-.0 

ei_l  =  D" 

en  _  -  w.  . 

i-i        l-i 


PMLi  - 1 


PJ:!?  =   1 
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STEr  5 


(DO     i-1,2     J 


Yes 


0^     =      Dn 


3UK3   -  SVK^   +  83 
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CONTINUE 


(j)0   i=l,  3^- 


<.  0 


S1M±  -  \'l± 


>0 


INFEASIBLE  REGION 


<  0 

| (  SUKW.  -  W 

l    l 


>  0 


INFEASIBLE  REGION 


kkk  =  0 


DO  i=l,4 


> 


Yes 


kkk  =  kkk+1 
mkckkk  =  n 


CONTINUE 
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Ye 

s 

ia  = 

=  i+1 

58 


59 


r 

v      YHS 

kc   =kc+] 
fflkkc   =  n 

CALL  AIL?l(kna,!anb,i,ia,T,D,TK) 


60 


III       =  kna 


CALIi  CAIII?(M, ,  ZOP, a, b,0,i ,  31,Hy) 


1 

^11  =  ^la1^ 
I.L       =  Jcab 

CAIL  CALliyCLLjZOP.a.b.e.i,  ;jl,Hy) 

YES 
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CALL  ALLT2(hma,kab,kmc,±,ia,Q,D,'.PK) 


- 

LL    =kna 

CALL  CALHV(LL,ZOP,p.(b,©,I,;jl,HT) 

HH7T,=HH7T-+HVVT 
Ij.L            JjX         Jj.Ij 

LL     =krab 

i 


CALL  CALHY(LL,ZOP,a,'b,©,i,Jl,H?') 


i 


LL     =kmc 


I 


CALL  CALI-IV (LL ,  ZOP ,  a,  b ,  6 ,  i  ,  jl ,  117) 


62 


HHVja  =  H^H+^n 


63 


.  -2bf  .«f-(af-af ) 


estJM  =  jf.-suKD, 

esui  =  V/.-3U19Q 
1    i 


YES 


Z  =  ZOP?1 


Z  .  _215iim.e3l.I.-(ailm-a^m) 


CALL  CALQ1  (Z,a,b,i,iin,ns,^IC) 


CALL  CALOl  (Z,a,b,i,im,ns,qiQ) 


im  =  0xm+1 
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CALL  CAL22 (a, b, i, in, jm, eSUM,ns,Z ) 
CAXI  CAI91(Z,a)b)i,im,ns,e>J9) 


o?m  =  wi-suMDi-e^m 


IX: 


Z  =  ZOPv1 


CALL  CAJZ2  (a, b,  i ,  in,  Jm, GSt'I  ,ns , Z  ) 

GALL  CAL9i(3,a,b,i,ira,ns,e,IG) 
GALL  CAiei(S,a,b,i,3n,nr,,e,I9) 


CALL  OAiei (Z , a , b , i , im, ns, 6, 19 ) 
CALL  CAL61  (Z  ,  a,  b ,  i ,  3m,HS ,  Q-,  10  ) 


G.:'m  =  GVm+1 
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esm.1  =  v/1-sro.H)£ 


jra  =  me2 
km  =  mo. 


CALL  GAlZ3(afb,ns,im,3m,km,©SUM,Z) 
I 


CALL  CAL91(Z,a,b,i,ira,ns,G,I9) 
CALL  CAL91(Z,a,b,i,  5m,ns,9,IG) 


CALL  CAI33(a,b,ns,im,  jm,km,eSUII,Z) 
CALL  CAL01 (Z , a, b, i, im,ns,©, I©) 

CALL  CAL91(Z,a,b)i,jm,ns,G,.I9) 
CALL  CALGl(Z,a,b,i,km,ns,G,IG) 


CALL  CAiei  (Z ,  a,  b ,  i ,  im,  ns ,  G ,  19 ) 
CALL  CALG1  (Z  ,  a,  b ,  i ,  jm, ns ,  G ,  16 ) 
CALL  CAL«l(Z,a,b,i,km,ns,9,I9) 


ei'-=er+i 


66 


67 


im  =  me. 

Jo  =  taep 

Ian  =  me. 

lm  =  ne. 

Zl 

^  •.  im  ,  im  ,km  .lm  „ 
=  2-bi   -bi   «bi   -bi   «Wi 

Z2 

^f.bf.bf.Caf-a^ 

Z5 

.  im  .km  .In  ,    in     jnu 
=  b±   .b±   .^   ,(aj  -a^  ) 

Z4 

=  bf  .bf  .bf  .(af-af ) 

Z5 

.   of  .bf  .bf^af-af ) 

•Z6 

.  jm  .km  .  lm,.  in  .1-an  .lm 
=  \   '\   '\  +°1   '\    -\ 

Z7 

=  bf  .bf  .bj 

B+bJ».b£".bfB 

Z  = 

■ 

(Z1+Z2+Z3+Z4 

+Z5)   /  (Z6+Z7) 

£m  -- 

.   (Z+af-af )    /   (2 

.b*»)     I 

e?rn  = 
i 

.   (Z+r 

.f-af )   /   (2.bf ) 

IMEGSRIZE 
> 

6^  = 

=   (Z+af-a^1)   /   (a.bf) 

&~    d:,:    £" 
wi    '    i    '    i 

#■' 

1111 

q  =q+1 


© 


BO.        /kfc=0\    Y  B 


3i"  =  W±— SDEIDj. 


© 

+ 


ime  =  me. 


e^e=-(Supf+ap-a^e)/(2.bf e)  ;  itoshize  ef 


69 


jme  =  me„ 


I 


9  fe=-  (lilpf +af  e-af  e )/  ( 2 .  bfe ) 
ef  e=-  (Eupf+af  e-af s  )/(2  .bf  e ) 


q  =  q+1 


eJn=\7i-su:.TD.-e^ne-0JI!le 


ime 

= 

He. 

3me 

= 

mep 

kme 

= 

me-. 

Gfe  =  -(ZUpin+a^!ae-a^me)/(2.t,|lae) 
efle  =  -(zWf+afe-aJme)/(2.b|me) 
63™e  =  -(ZUpin+afle-afle)/(2.brne) 

©I™    =  \Yi-e^me-Q^m8-Qtjne 


HWSSERIZE 

ekme 


20] 

f 


a  =  ot1 
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21 

T 


irae  =  ne 


ke 


ef1G 


- (ZUP^n-l-a^UL-a^ie)/(2 .b^1 1 " ) ; IMEGEHIZE  0 


ime 

i 


1 

ime 

= 

me1 

jme 

= 

iae2 

iwe  _ 
wi       - 

-(ZU?f-iafe-afe)/(2.bfe): 

UJ33GSRIZE 

of e  = 

-  (znpf1+4ae-a3me  )/(2 .  bf  e ) 

pirae-jne 
Hi      »*i 

q   =   q-i-1 
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72 


73 


74 


SU1'95=0 


("lO  k  =    1,2\*- 


5*  =  Dn-9*-G* 


75 


Step  7 


SUIlF  =  0 


r 


(j>0  it  ---    1,?} 


1 

[TO   i  =   1,?)- 


* 

I?  .  a^^Ce?)^ 

suiicn  =  sieic'W.1 

COETIHUS 


SUIiP  =  SUKP+SDKC 


OOMHTOE 


GO   TO  SI3E  8 


STEP  8  &  9 


SHIN  =  SUT-'.F 


jomin  =  jo 
jmin  =  j 


f  DO  n 

-      1           '^-1 

1 

PUNCH  e£,  ©g,   «" 

76 


CONTINUE  V 


PUNCH  jomin,  jmin,  SKIN,  SUKF 


GO  BACK  TO  STEP  5  FOR  A  NEW  REGION  OF  SEARCH 
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SUBHOUTIiTE  ALIS1  (kma,tanb,i,ia,e,D,0!H) 

DIHISSIOH  e(3,4),D(4)!a,H(2), 


_kma       _ima  ,,Iaaa 
9i        ~  D       ~faia 


110 


Y3S 


s.    u           xa 

ef *  . 

6f  <  % 

ekmb_  0 

a. 

HETDHN 

SUBROUTINE  CALHV(LL,ZOP, a,  15,6,1,  j,HV) 
DBCWSIOH  20P(2,8),a(3,4),T5(3,4),ei3,4),Hv(3) 

irvL1  =  (zopl+a^'-a^.el^+b^cejp')2 

B2OTEH 

EI'ID 

7S 


SUBROUTINE  ALLT2  (fena,  tab,  tec ,  j  §  la>  e>  D>  ^ 
DIMENSION   0(3,  2j.),  D(4)t  TH(2) 


e^3-  =  Dkma  -  ekma 


ia 


YES 


1%   =  S%    -   G^113- 
3kmb  b  Dkmb  _     kmb 


ekmb  =  0 

ef 10  .  o 


YES 


,fcoc 


=  ra4  -  e*mb 


RETURN 
END 
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SUBROUTINE. 
DIMENSION 

CAIOl  (Z,a,bfi,im,ns,e,ie) 
a(3.4)i'bC3,4),e<3,4,), 10(3,4) 

1 

1 

ief= 

-(Z+aJ"«g)  /  2.bf 
-(Z+aJ>g)  /  2.bf 

1 

EO 

0ira 
l 

= 

IOim 
1 

SUBROUHHB     CAXZ2    (a,b,i,iin, jia,OSUI-I,ns,Z) 
DIKEKSIOH       a(3,4),b(3,4) 


1 


Zl  =  -2.bf\b?'m.e3UK  -  b?m(afm  -  aim) 
11  1       i  ns' 


Z2  .  bj* (af  -  aJ°) 


Z  =   (Zl  -   Z2)/(b^m+  bj>m) 


RBTOHH 
END 


&0 


SUBBOraHTB  GAI23  (a,b,ns,  ira,  jm,im,eSDH,Z ) 
DUEBSIOH  a(3,4),   b(3,4) 


zi  =  2bf.bf.bf.esuH 

Z2  -  bf  ,bf  (af-af ) 
Z3  =  bf  .bf  (af-a-g) 
24  =  bf  .bf  (af-a^) 
35  =  bf  .bf +bf  .bf +bf  .bf 


7   -  -(Z1+Z2-:-33+Z4) 

z5~ 


R3TOKI 
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SUBROUTINE  AM3  (6 ,  i,tf  ,3D,  to,  la) 
DIMENSION  e(3,4),W{3):,D(4) 

T 


Ta=Tir-i 


NO 


YSS 


roiai 


END 
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APPENDIX  B  SYMBOL  TABLE  FOR  COMPUTER  PROGRAM 


Program 
Symbol 


Mathematical 
Symbol 


Explanation 


A(I,N) 
B(I,N) 
BB(I,J) 


CKK 

CO 

COUNT 

D(N) 

F(I,N) 


F^(e^) 


H(I),HV(LL),HV(L1)   H^ 
I  i 

IA 


The  coefficient  of  the  first 
order  theta  in  the  cost  func- 
tion expression. 

The  coefficient  of  the  second 
order  theta  in  the  cost  func- 
tion expression. 

Same  as  above  but  the  super- 
script is  ordered  according  to 
the  value  of  the  upper  break- 
ing value. 

A  counter  to  indicate  whether 
the  8j  =  equation  of  optimality 
have  been  allocated.   It  is  use 
for  allocating  the  severely 
constrained  condition. 

A  counter  of  the  number  of  the 
states  with  the  severely  con- 
strained condition. 

A  counter  to  indicate  whether 
the  ZOPJ°  (3c)  is  out  of  order. 

The  number  of  time  step  6  is 
repeating  when  there  is  a 
severely  constrained  condition. 

The  number  of  units  of  the  re- 
source required  by  the  rr1"  de- 
mand point. 

The  cost  incurred  by  transport- 
ing e£ 

The  Hamiltonian  function. 

Subscript  for  origin.  ■ 

Subscript  of  the  next  or  pre- 
ceding origin. 
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IM,  IKE 

IMK 

IT 
J 

JM,  JME 
Jl 

jmin 

JOMIN 
K 

KC 

KE 

KK 
KKK 

KMA 


The  index  of  the  first  stage 
to  have  a  6"  =  equation  of 
optimality  condition. 

The  index  of  the  first  stage 
to  have  a  severely  constrained 
condition. 

The  integerized  9^. 

The  subscript  for  the  .ordered 
region  of  search,  ZOPJ. 

The  index  of  the  second  stage 
to  have  a  6"  =  equation  of 
optimality  condition. 

Superscript  of  the  ordered  re- 
gion of  search,  ZOPJ1. 

The  optimal  region  of  search . 
of  the  second  origin. 

The  optimal  region  of  search 
of  the  first  origin. 

A  counter  for  the  number  of 
free  variation  of  6n  condition; 
It  is  also  used  for  subscript 
in  place  of  i  for  a  loop  in- 
side i-"DO-loop". 

A  counter  for  the  number  of 
free  variation  of  6£  condi- 
tion. 

A  counter  for  the  number  of  6?  • 
equation  of  optimality  con- 
dition. 

A  counter  for  the  number  of 
severely  constrained  condition. 

A  counter  for  the  number  of 
stages  in  the  origin  which  has 
Dn  larger  than  W±. 

The  index  of  the  first  stage 
to  have  a  severely  constrained 
condition. 


fit 

p  A  counter  to  indicate  whether 

the  61  with  the  severely  con- 
dition has  been  allocated. 

PM(I,N)  Identification  for  stages  in 

the  origin  where  Dn  is  larger 
than  Wi. 

Q  A  counter  to  indicate  whether 

the  9j  with  the  free  variation 
of  6n  conditon  has  been  allo- 
cated. 

SMIK  Cgjj       The  current  minimum  total  cost 

of  transportation. 

SUM(I)  x£        The  state  variable,  the  total 

amount  of  resources  which  have 
been  transported  from  the  i*" 
depot  to  the  first  n  stages 
(demand  points). 

SUM  C(R)        £  F?(©J)       The  total  cost  of  transpor- 

i=l  ting  6"  from  all  depots  to  the 

n   demand  point. 

SUMD(I)  Similar  to  SUM(I)  exceDt  only 

8n  =  Dn  are  added. 

SUM  F  %s  The  total  cost  of  transporting 

0j  from  all  depots  (origins) 
to  the  first  n  demand  points. 

SUM  TT(I)  Similar  to  SUM(I )  except  only 

0n  =  Dn  and  60?  =  Dn-0ja  (where 

6ia  =  equation  of  optimality 
and  06la  has  a  severely  con- 
strained condition). 

SUM  T(3)  Similar  to  SUM(I)  but  1-3 . 

SUM  W(I)  The  constraint  for  the  stages 

with  0n  =  0. 

TU>N)  %        The  quantity  of  the  resources 

sent  from  the  ith  depot  (origin) 
to  the  ntfl  demand  point. 

TH(I)  The  available  resource  to  be 

allocated  to  the  free  varia- 
tion of  9j  condition. 
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TR 
18(1) 

TSU 
TSUM 

TT(I,N) 

TM1,  TM2,  TM3 

W(I) 


w. 


3      U 
WT         Z  wi  =  E  Dn 
i=l     n=l 


ZUP(I,N),ZBU(I,J) 

ZOP(I,J) 
ZZ 


Zl,  Z2,  Z3,  Z4,  Z5, 
Z6,  Z7 


7n 


2L0W(I,N),ZBL(I,J)   z" 


The  remainder  of  9i  in  the  origin 
due  to  integerization  of  9?. 

The  available  9?  in  the  origin 
when  maximum  available  Gi'.is 
allocated  to  the  severely 
constrained  condition. 

Same  as  TS(I). 

The  resource  available  for  8t   = 
equation  of  optimality  condition 
in  the  state. 

The  9^  in  case  where  maximum 
available  9?  is  allocated  to 
the  severely  constrained 
condition. 

The  9i  with  the  severely  con- 
strained conditions  which  give 
the  minimum  value  for  the  Ham- 
iltonion  function. 

The  total  number  of  units  of 
the  resource  available  at  the 
ith  depot. 

The  total  resource  transported. 


The  value  of  zf  used  in  finding 
9.  =  equation  of  optimality. 

The  lower  limit  of  the  breaking 
value,  Zj,  which  is  obtained 
by  equating  the  derivative  of 
Hy  with  respect  to  9n. 

The  upper  limit  of  the  breaking 
value,  Z§\ 

The  regions  of  search. 

The  dummy  of  Zj°  value  from  the 
previous  loop  of  Step  6. 


The  terms  of  the  equation  for 

calculation  of  the  Z"  value. 

i 
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KMB  '  The  index  of  the  second  stage 

to  have  a  severely  constrained 
condition. 

KMC  The  index  of  the  third  stage 

to  have  a  severely  constrained 
condition. 

KME  The  index  of  the  fourth  stage 

to  have  a  @P  =  equation  of 
optimality  condition. 

L  i         Index  for  i. 

LL,L1  Indexes  for  stages  with  severely, 

constrained  condition. 

M(K)  Indexes  for  stages  with  free 

variation  of  0n  condition. 

ME(KE)  Indexes  for  stages  with  9n  = 

equation  of  optimality  con- 
dition. 

MK(KC)  Indexes  for  stages  with  severely 

constrained  condition. 

N,  Nl  n        Superscript  of  the  demand  points 

or  stages  n  =  1,  2,  ...,  N. 

NJ  Total  number  of  regions  of 

search  for  the  second  origin. 

NJO  Total  numbers  of  regions  of 

search  for  the  first  origin. 

NK(I)  Total  number  of  regions  of  search 

in  each  origin. 

NN  N        Total  number  of  the  demand  points 

or  stages. 

NO  The  index  of  the  largest  break- 

ing value  in  the  origin. 

NS 


Total  number  of  the  depots. 

NJ1  Total  number  of  regions  of 

search  in  each  origin. 
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APPENDIX   C    COMPUTER  PROGRAM 


//MASTER        JOB   00I50406G0C1 , SKONGKATONGS  20,20,900 

0  0 

//    EXEC    FTECLGKS 

//FORT.SYSLIN  DO  UNIT=194 

//FORT. SYS  IN    CD    » 

DIMENSION  ZGP(2,8),ZUI3,4),ZBL(2,4),NK(2),W(3),D(4) 
DIMENSION  H ( 2 ) , SUM ( 3 ) , SUMW ( 3 ) , TT ( 3 , A ) , SUMO ( 2 ) , ME ( 4 ) 
DIMENSION  A(3,4),P(3,4) ,T( 3,4) ,BB 1 3,4 ) , ZUP ( 3,4) 
DIMENSION  SUMTT(2),PM(  3,4  )  ,MPC  (4)  ,  Z10M  3,4  ) ,MK(4) 
DIMENSION  TH(2),SUMC(4),F(3,4),IT(3,4),SUMT(3),TS(2) 
DIMENSION  HHV(6),ZZ(2),M(4),HV(3),ZBU(2,4) 

100  F0RMAT(3F8.4,3F10.0) 

102  FURMAT(2HI=,I3,5H    J=,I3,16H        ZBU I I , J )=,F9. 5 ) 

103  FORMAT! 28HTHIS  IS  AN  INFEASIBLE  REG] 

107  FORMAT! 14) 

108  F0RMATI3F15.3) 

101  FQRMAT(2HI=, I3.4H   N=,I3,11H   ZUP  (  I  ,  )-,'   .3, 
L12H   ZL0W(I,N)=,F9.3) 

104  FORMAT(3HJ0=,I2,5H    J=,I2,15H    ZOP ( 1-1 , JO )=,F9.3 , 
C12H    ZOP! I,J)=,F9.3) 

109  F0RMATI8H   J0MIN=,I2,8H    JMIN=,I2,9H     SMIN=,F10.2, 
C9H     SUMF=,F16.2) 

110  FORMAT! IX, 11HTH1RDCHECK=, 2F9.2 ) 

111  FORMATIF10.3) 

112  F0RMATI2F15.6) 

113  F0RMATI2F9.1) 

114  F0RMATI1X,11HFIRSTCHECK=,4I3,9F6.1 ) 
L0=0 

NN=4 
NS  =  3 
900  DU  1  1=1,2 
DO  2  N=1,NN 
READ!!  ,100)  A  (  I  ,  )  ,  (I  ,  ,  )  ,  A(  NS  ,  N)  ,  D(  N)  ,H  I  I )  ,  W  (  NS  ) 


IF(B(I,N))40,5C, 

40  T(I,N)=D(N)/2. 
ZUPU,N)=-A<I,N)+A(NS,N)-2.»8< I.N)*T(I.N) 
ZL0W(I,N)=ZUP( I.N) 

GO  TC  2 

50  ZUPUfN)=-A(I,N)+A(NS,N) 
ZLOW(I,N)=ZUP(I,N) 

GO  TO  2 

60  ZUP<  I,N)  =  -AU,N)+A(NS,N) 
ZLOW(I,N)=-A(I,N)  +  A(NS,M-?.*R(I,N)»D(N) 

2  WRITE(2,101>I,N,ZUP(I,N),ZL0W(I,N) 
SMIN=1000000. 

DO  3  J=1,NN 

Z8U(I,J)=-1000. 

00  4  N=1,NN 

IF(ZBU(I,J)-ZUP(I,N>)41,51,61 

41  ZBU(I.J)=ZUPII,N) 

51  NO=N 

GO  TL  4 

61  GO  TO  4 
4  CONTINUE 

ZU(I.N0)=ZUP(1.N0) 
ZUP( I,NO)=-1000 
B8(I,J)=B( I, NO) 
ZBLU,J)  =  ZLOW(I,NO) 

3  HRITE(2,lQ2II,J,ZBUtI,  J) 
K=o 

00  5  J=1,NN 
IF(BR< I, J) 142,52,62 

42  K=K+1 
ZOP(I,K)=ZBU(I,J)-.Q 
GO  T 

52  K=K+1 
ZOP(I,K)=ZBU(I, J) 
K=K+1 


89 


ZOP(I,K)=ZBU( I, J)-. 005 
GO    TC    5 
62    K=K  +  1 

ZOP(I,K)=ZBU( I, J)-. 005 
K=K+1 

ZOPU,K)  =  ZBL< I, J)-. 005 
5    CONTINUE 

NK(I)=K 
1    WFUTE(2,1C7)NMI> 
NJO=NKI  1) 
NJ=NK(2) 
1  =  1 
NJO=NJO-1 

316  C0=0. 

00  300    J0=1,NJ0 

IFIZ0PI I,J0)-Z0P( I.J0  + 1)1340,360,360 
340    ZAP=Z0P(I,J0) 

ZCP(I,JC)  =  Z0PU,J0  +  1) 
ZuP( I,J0+1)=ZAP 
C0=C0+1. 

360  GO   TO    300 
300    CONTINUE 

IF(C0)361,351,361 

361  GO  TO  316 
351  NJ0=NJ0+1 

DO  500  JO=1,NJO 
500  HRITE(2,111)Z0P(I,J0> 

1  =  2 
NJ=NJ-l 

317  C0=0. 

DO  301  J0=1,NJ 

IF(Z0P(I,J0)-ZOP(  I,JQ+1>)  343,363,363 
343  ZAP=Z0P(  1,J0) 

zupii, jo)=zop(i, ju+i) 

ZGP1 I.JL+1 )=ZAP 
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C0=CC+1. 

363  GO  TC  301 
301  CONTINUE 

IF(C0)364,354,364 

364  GO  TC  317 
354  NJ=NJ+1 

DO  501  J0=1,NJ 
501  WRITE(2,111)Z0P(1,J0) 
WT=W(1)+W(2)+W<3) 
DO  6  J0=1,NJ0 
DO  8  J=1,NJ 
L=l 
I  =  L+1 
Q=0. 
C=0. 
P=0. 

COUNT=  . 
CKK=C. 

ZZI21-  . 
DO  7  K=l,3 
SUM[K)=0. 
SUHW(K)=WT 
DO  800  N=1,NN 
800  PM<K,N)=0. 
7  CONTINUE 

Wr<ITE(2t104)J0,J,ZGP(I-l,J0>,Z0P(I,J) 

DO  9  N=1,NN 

ZUP( I,N)=ZU< I,.) 

ZUP<I-1,N)=ZU( I-l.N) 

IFIZOPI I, J 1-ZUPI1, v 1)4  5,55,65 

45  IF(ZOP(I,J)-ZLOWU,N)  )  46,  56,  56 

46  IFIZCPI  I-l,JO)-ZUP(  I-l.N)  >47,r>7,57 

47  IFIZCPI I-1,JO)-ZLOHI I -  1 ,N ) ) 48, 58, 58 

48  H(I-l)=(ZGP(  l-l.JOl+AI  I-1,N)-A(NS,N))»U(N)  +  I  1  - !  ,  J )»D I N)«»2 
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H(I)=(Z0P(I,J)+A(I,N)-A(NS,N))*D(N)+B(I,N)»0(N)««2 
IF(H(I-1)-H( [ 1)49,59,59 
49  IF(D(N)-W(I-1))830, 830,850 

830  T(  I,N)=0. 
T(  I-l,N)  =  l)(N) 
GO  TO  116 

850  T( I,N)=D(N)-W(I-1) 
T(I-1,N)=W(I-1) 

801  PM(I-1,N)=1. 
PM( I,N)=1. 
GO  TC  116 
59  IF(D(N)-W(I))831«831,851 

831  T(I-1,N)=0. 
T(I,N)=0«N) 
GO  TC  116 

851  T(I,N)=W( 1) 

T( I-1,N)=D(N)-W( I) 

GO  TO  801 
58  T( I-1,N)=-(Z0P< I-1,J0)+A(I-1,N)-A(NS,N) )/ ( 2 .*B ( 1-1 ,N) ) 

T( I,N)=D(N)+2. 

GO  TC  9 
57  IF(D(N)-W(I))832,832,852 

832  T(I,N)=D(N) 
T( I-1,N)=0. 
GO  TO  116 

852  T(I,N)=W(I  ) 
T(I-1,N)=0. 
GO  TO  801 

56  IF(ZOPII-l,JO)-ZUP(l-l,N))70,80,90 

70  IFIZOPI I-1,J0)-ZL0H( I-1,N))71,81,81 

71  T( I,N)=-(ZOP( I,J)+A(I,N)-A(NS,N))/(2.*B(I,N)) 
Tl  I-1,N)=D(N)*2. 

GO  T 
81  T(I,i\)  =  -(ZOP(I,J)  +  A(I,N)-AINS,N))/(  .  ■>  (  I  ,  N  )  I 

T ( 1-1, N )=-( ZOP( 1-1, JO ) +A( I-1,N)-A(NS,N) 1/(2. •8(1-1, N)) 
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GO  TO  9 
80  T( I,N)=-(ZOP(I,J)+AU,N)-A(NS,N))/(2.»B(I,Nl) 

T< I-1,N)=D(N)+1. 

GO  TO  9 
90  T(I,N)=-(Z0P(I,J)+A(I,N)-A(NS,N))/(2.«B(I,N)1 

T(  1-1,  N)  =  .  . 

GO  TO  116 
55  IF(ZaP(l-l,J0)-ZUP(I-l,N))72,82,92 

72  IF ( Z0P( 1-1, JO )-ZLOW( I- 1,N 1)73,83,83 

73  1F(0(N)-W( 1-11)833,833,853 
833  T(I,N)=0. 

T( 1-1,N)=D(N) 
GO  TO  116 
853  T(I,N)=0. 

T(I-1,N)=W(I-1) 
GO  TO  801 

83  T(I,N)=D(N)+1. 
TU-1,N)=-(ZOP(I-1,JO)+A(I-1,N)-A<NS,N)  )/< 2 .»Bt I-l ,N) ) 
GO  TO  9 

82  T(I,N)=D(N)+1. 

TII-1,N>=D(N)+1. 

GO  TO  9 
"2    T(I,N)=D(N)+1. 

T( I-1,N)=C. 

GO  TO  116 
65  IF(ZCP(I-1,J0)-ZUP(  I- 1,iM  1)74,84,94 

74  IF(ZCP( I-l.JO  WLOWJ I- 1,N) 175,85,85 

75  GO  TC  73 
85  T( I,N)=0. 

T( I-1,N)=-(Z0P< 1-1, JO) +A(I-1,N)-A( NS.rO  )/ ( 2.*3( 1-1 ,N) ) 
GO  TC  116 

84  T(I,N)=0. 
T(I-1,N)=D(N)+1. 
GU  TC  116 

94  T(1,N)=0. 
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T(1-1,N)=0. 
T(NS,N)=D(N) 
GO  TO  117 

116  DO  10  K=l,2 

I  F  (  PM  (  K  ,  N  )  - 1 .  )  8  34 ,  844  ,  8 34 

834  IFIT(K,N)-0(N))76,86,76 
844  IF(T(K,N)-W(K))835,86,835 

835  IF(T(K,N))836,87,836 

836  IF(T(K,N)-D(N))846,10, 10 
846  SUM(K)=SUM(K)+T(K,N) 

GO  TO  10 

76  IF(T(K,N))77,87,77 

77  GU  TO  10 

87  SUMW(K)=SUMW(K)-D(N) 

GO  TO  10 
86  SUM(K)=SUM(K)+T(K,N) 

SUMW(3)=SUMW(3)-T(KfN) 
10  CONTINUE 

GU  TO  9 

117  T(3,N)=D(N) 

SUM(3)  =  SUM3)+T(3,M 
9  WRITE(2,112)T(1,N),T(2,N) 
00  11  K=l,3 
IF(SUM(K)-W(K))78,78,98 

78  IFISUMWIK)-W(K) )79,   , 

79  WRITEI2.103) 
IX=1 

DO  405  N=1,NN 
WRITE(2,113)T(1,N),T(?,N> 
405  CONTINUE 

WRITE(2,107)IX 
GO  Ti 
89  KKK=0. 

DO  802  N1=1,NN 
IF(T(K,N1> 1838,802,838 
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838  KKK=KKK+1 

MPC(KKK)=N1 
802  CONTINUE 

IF(KKK-1)11,849, 11 
849  N=MPC(1) 

IFIPIMK,  N)  1840,11,840 
840    IF(K-1)870,880, 

870  KA=K-1 

GO  TC  810 
880  KA=K+1 
810  IF(T(K,N)-H(K))871,11,  871 

871  IF(D(N)-W(K))11IU,882 
882  T(K,N)=W(K) 

T(KA,N)=D(N1-W(K) 

GO  TO  11 
98  GO  TO  79 
11  CONTINUE 
200  DO  12  1=1,2 

SUM(I}=0. 

SUMD(I)=0. 

SUHTT(I)=0. 

K=0 

KE=0 

KK=0 

00  13  N=1,NN 

IF(PM  I.N1-1.  1837,847,837 
837    IFITI I, N)-< DIN  1  +  1.  1)1 V   ,150,160 
847    GO    TO    151 

140  IF(T(I,N)-(D(N) 11141, 151, 151 

141  IF(T(I,N))142,152,  162 

142  GO    TC    79 
152    GO    TO    118 
162    KE=KE+1 

HE(KE)=N 
118    SUM! I)=SUM( Il+T ( I ,     ] 
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GO  TC  13 
151  SUMD(I)  =  SUMO( I)+T(  1,M) 

SUMTTI I)=SUMTT(  I)  +  T( I,   ) 

GO  TC  118 
150  K=K+1 

M(K)=N 

GO  TC  13 
160  KK=KK+1 

IF(I-1)143,153,143 

143  IA=I-1 
J1  =  J 
NJ1=NJ 

119  IF(C) 144,154, 144 

144  IF(TS(I))347,357, 
347  GO  TC  13 

357  T(I,N)  =  D(N)-TUA,N) 

P=P-1. 

GO  TC  151 
154  TT(  I,N)  =  D(N)-T( IA,N) 

P=P+1. 

SUMTTI  I)  =  SUMTTt  Il+TTI I,N) 

GO  TC  13 
153  IA=I+1 

Ji  =  JC 

NJ1=NJ0 

GO  TC  119 
13  CUNTINUE 

IFIC1145, 155,145 
145  IF(KK)34S,348,368 

348  GO  TO  I? 

368  IF(TS(I))349,359,359 
359  Gu  TC  12 

349  KC=0 

TH( I)  =  w( I) -SUM! 1) 
IFITHI I  1)370,38  , 
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370  GU  TC  79 

380  DU  303  N=l,4 

IF (T ( I, N)-( DIN  1+2. 1)371,381,371 

371  GU  TL 

381  KC=KC+1 
MK(KC)=N 

303  CONTINUE 

IF(KC-2)373,3R3,393 

373  IMK=WK(1) 

T(  I,  IMK)=TH( 1 1 

P=P-1. 

GO  TO  12 

P=P-2. 

DU  304  Ll=l,2 

IFILl-21374,384,374 

374  KMA=MK(11 
KHB=«K(2) 
GO  TO  319 

384  KMA=KK<2) 

KMB=MK( 1) 
319  WRITE(3,114)KMA,KMB,I, IA,T<  I, 1 ),T(  [,,)  ,  I  (1  ,  ),r(l,4), 
CD(1),D(2),D(3),IJ(4),,TH(  J) 

CALL  ALLT1  ( KMA.Kf K, I , I A,T, 0, TH) 

WRITE(3,11C)T(I,KMA),T(  I.KHB) 

LL=KMA 

HHV(L1)=0. 

CALL  CALHV  (LL , 20P , A, B , T, I , J  1 ,HV 1 

HHV(L1)=HHV(L1)+HVILL) 

LL=KNB 

CALL  CALHV  ( LL, ZOP , A, H , T, I , J  1 , HV 1 

HHV(LI)=HHV(L1)+HV(LL) 

IFUl-21375,385,375 

375  TM1=T(I,KMAJ 
TM2=T( I,KMB) 
GO    TC    304 
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304    CONTINUE 

385  1FIHHVI 1)-HHV(2))376,3S6,386 
J76    T(I,KHA)=TH1 

T(  I,KMB)=TM2 
GO    TO    12 

386  GO    TO    12 
393    P=P-3. 

DO    305    Ll=l,6 
IF(Ll-2)377,387,397 

397  IF(Ll-4)378, 388,398 

398  IFUl-61379,389,389 
389    KMA=MK(3) 

KMB=MK<2) 

KMC=VK(1) 

GO  TO  320 
379  KMA=MK(3> 

KMfi=MK(l) 

KMC=MK(2) 

GO  TO  320 
388  KMA=MK(2) 

KMB=«K(3) 

KMC=MK( 1) 

GO  TO  320 
378  KMA=MK(2) 

KMB=MK(1) 

KMC=MK<3) 

GO  TO  320 
387  KMA=MK(1) 

KMB=MK(3) 

KMC=MK(2) 

GO  TO  320 

377  KM4=*K(1) 

-MKI2) 

KMC=MK(3) 

GO  TO  320 
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320  CALL  ALLT2  ( KMA ,KMB,KMC , I, I A,T, D, TH) 
HHV(L1)=0. 
LL=K*A 

CALL  CALHV  ( LL , ZOP , A, B , T, I , Jl , HV ) 
HHV(L1)=HHVIL1)+HV(LL) 
LL=KMB 

CALL  CALHV  ( LL , ZQP, A, B , T, I , J  I, HV ) 
HHV(L1)=HHV(L1)+HV(LL) 
LL=KMC 

CALL  CALHV  ( LL, ZOP , A, B , T, I , J  1 , HV) 
HHV(L1)=HHV(L1)+HV(LL) 
IF(L1-2)440,45C,V 

440  TM1=T(T,K«A) 
TM2=T(I,KMB) 
TM3=T(I,KMC> 
GO  TO  305 

450  IF(HHV(L1-1)-HHV(L 1)1441,451,451 

441  T(I,KMA)=TM1 
T(I,KMB)=TM2 
T(  I,KHC)  =  TM3 

451  GO  TC  305 
305  CONTINUE 

GO  TC  12 
155  IF(K-1)147,157,167 
147  IF(KE-1)148,158,168 
14B  GO  TC  1? 
158  IM=ME(1) 

IF(KK)470,470,490 
470  T( I,IH)-H( I)-SUMD( I) 

Z=-2.»R(I,IM)»T(I, IH)-(A(I,IM)-A(N5,IM) ) 

744  IFIZ-ZCPI I,J1))745,79,79 

745  IF( J1-NJ1)746,12,12 

746  IF(Z-Z0PU,J1  +  1))79,79,12 
490  TSUM=HU)-SUMD<  I) 

TSUI=W{1)-SUMTT(I) 
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IF(TSUI)471,481,481 

471  Z=Z0P(I,J1) 

CALL  CALT1  (Z,A,B, t,KM,NS,T, IT ) 

TU,IM)=T< I,IM)+1. 

GO  TO  323 
481  Z=-2.«8(I,IM)«TSUI-(A< I , IM )-A(NS> I M) ) 

CALL  CALT1  < Z , A , B, I , IM, NS, T, IT ) 

GO  TC  323 
168  IF(KE-3)149,159,169 
L4)  TSUM=W(I)-SUMD(I) 

IM=ME(1> 

JM=ME(?) 

IF(KK)472,472,492 

472  CALL  CALZ2  ( A, B, I , IM, JM.TSUM.NS , Z ) 
CALL  CALT1  ( Z , A ,B, I , IM, NS, T, I T ) 

T(  I,  Jl*)=W< I)-SUMD( I)-T(  I.IM) 
GO  TO  744 
492  IF(TSUI>473,483,483 

473  Z=Z0P(I,J1) 

CALL  CALT1  t Z , A,B, I , I M , NS, T, I T ) 

CALL  CALT1  ( Z , A,B, I, JM.NS.T, IT ) 

T( I, IM)=T( I,IM)+1. 

T(I,JK)=T{ If JM)+1. 

GO  TO  323 
483  CALL  CALZ2  ( A,B , I , IM, JM.TSUI ,NS, Z ) 

CALL  CALT1  ( Z , A, B, I, IM , NS, T, IT ) 

CALL  CALT1  ( Z , A ,B, I , JM , NS,T, I T ) 

GO  TO  323 
159  TSUM=W( I)-SUMD( I) 

IM=ME(1) 

JM=ME<2) 

KM»M  (  ) 

IF(KK)474,474,494 

474  CALL  CALZ3  I  A, B ,NS , IM, JM.KM, TSUM, Z  ) 
CALL  CALT1  (  Z  ,  A  ,B  ,  I ,  I  M,  NS,  T,  I  T  ) 
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CALL  CALT1  ( Z, A.B, I, JM, NS, T, I T ) 

T(I,KM)=W< I)-SU«D( I)-TU,IM)-T(I,JM) 
'   744 
494  IF(TSUI)475,485,485 
475  Z=Z0P(I,JI) 

CALL  CALT1  (Z, A,B, I, IM ,NS,T, IT) 

CALL  C'LTl  (Z,A,B, I,JM,NS,T, IT) 

CALL  CALT1  ( Z, A ,B, I ,KM , NS, T, IT ) 

TCI,IM)=T< I,IM)+1. 

T( l,JM)=T(I,JM)+l. 

T(  I,KM)=T<  I,KK)+1. 

GO  TO  323 
485  CALL  CALZ3  ( A,B,NS, IM, JM.KM, TSUM,Z ) 

CALL  CALT1  ( Z , A.B, I , IM , NS, T, IT ) 

CALL  CALT1  ( Z , A ,B, I , JM , NS, T, I T ) 

CALL  CALTI  ( Z, A, B, I, KM, NS, T, I T ) 

GO  TO  323 
323  IF(CCUNT)445,445,465 
445  CKK=CKK+1. 
465  IFIZZ1I )-Z)447,457,447 

447  ZZ(I)=Z 
GO  TC  12 

457  IF(Z-Z0P(I,J1))448, 448,468 

448  IF(J1-NJ1)543,469,469 

543  IF(Z-ZOP( I, J  1  +  1 1)468,468,469 
469  TS(I)=TSUI 

CKK=CKK-1. 

GO  TO  12 
468  GO  TC  79 
169  IM=ML(1) 

JM=ME(2) 

KM=ME(3) 

LM=ME(4) 

Z1=2.»B(I,IM)*B(I,JM)«B(I,K*I)«B(I,LM)«WII) 

Z2=B(I,JM)*B(  I,KM)»iJ(  I,LM)»(A(  I , IM )-A (NS, I M ) ) 
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Z3=B t I* IM)*B« I,KM)*B( I , LM) » ( A ( I , JM )-A (NS, JM ) ) 
Z4=B<I,IM)«B( I,JM)«B( l,LM)»(A(I,KM)-A<NSfKM>) 
Z5=B(I,IM)*B( I, JM)»B( I,KM)»(A(I,LM)-A(NS,LM1) 
Z6=R(I,JM)*B( I,KM)»B( I,LM)+B( ],IM)»B< I , KM) «3 ( I , LM) 
Z7=B(I, IM)»B( I, JM)«B( I,LM)+B( I,IM)*B(I, JM)*B(I,KM) 
Z  =  -Ul  +  7  2  +  Z3  +  Z4+Z5)/(Z6  +  Z7) 

IT(I,IM)=-(Z+AI 1,IM)-A(NS, IM))/(2.»P( I, IM) ) 
T(I»IHWT(Ii  IM) 

IT(I,JM)=-(Z+A(I,JM)-A(NS,JM))/(2.»B(I,JM) ) 
TCIf JM)=IT(I,JM) 

IT(I,KM)=-(Z+A( I,KM)-A(NS,KM))/(2.»BU,KM) ) 
T(I,KM)=IT(I,KM) 

T(  I,LM)=W( Il-TI I, IM)-T( I,  JM)-T(  I  ,KM) 
GO  TO  12 
157  IM=M(K) 

IF(KE-1)170,180,190 

170  IF(KK)171,181,171 

171  Q=Q+1. 
GO  TO  12 

181  T(I.IM)=W(I)-SUMD(  I) 
GO  TC  12 

180  IME=ME(1) 

IT(I,IHE)=-(ZUP(I, 1M)+A(I,IME)-A(NS,IME))/(2.»B(I,IME)1 

T(I,IM£)=IT(I,IME) 

IF(KK)172,182,172 

172  0=0+1. 
GO  TG  12 

182  T(It IM)-M( I)-SUMD(  I)-T(  I, IME) 
GO  TO  12 

190  IF(KE-3)173,183,183 

173  IME=ME(1) 
JMS=ME(2) 

IT(I,IME)=-(ZUP(I,IM)+A(I,IME)-A(NS,IME))/(2.»B(I,I«EJ1 
T(  I,IME)=IT< I, IME) 

IT(I,JME}=-<ZUP(I,  IMJ  +  ACI,  JME)-A(NS,JME))/(  •'.*  M  I,  JME)  > 
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T(I,JME)=IT( I.JMEJ 

IF(KK) 174.184.17A 
174  0=0+1. 

GO  TO  12 
184  T<I,IM)  =  W<I)-SUMD( I )-T (  I, I  ME >-T<  I ,  JME  ) 

GO  TO  12 
183  IME=ME(1) 

JME=PE(2] 

KME=ME(3) 

II (I,IME)=-(ZUPU,IM)+A(I,IME)-A<NS,IME))/{2.»B(I,IME) ) 

T(  I,IM  )-IT(I,IME) 

IT( I, JME)=-(ZUP(If IM)+A(I, JME)-A(NS,JME))/(2.»B(I,JME> ) 

T(I, JWt)=ITII,JME> 

IT(I,KKE)=-<ZUP(I,IM)  +  MI,KME)-A(NS,KMEn/(2.»R(I,KME)) 

T(  1,KME)=IT( I.KME) 

T( I ,IN)  =  W( II-TI  I,1MC)-T(I, JHEJ-TCI ,KNE) 

GO  TO  12 

IF(KE-1)175,185,195 

0=0+1. 

GO  TO  12 

IME=FE(U 

1T(I,IME)=-(ZUPII, IM)+A(I,IM£)-A(NS,IME))/I2.»B(I, I   ) ) 

TCI.IHEt-ITC I.IME) 

GO  TO  175 

IME=ME(1) 

JME=ME(2) 

IT(I,IME)  =  -(Zl)P(l,IM)  +  AU, IME)-A(NS,IME)>/(2.»8(I,IME)) 

T ( I , I K  f )  =  I T I  I.IHE) 

IT(I,JME)=-(ZUP(I,  I  M)  +  ACI,  JME) -AINS,  JME  ))/(.'.*  M  I,  JME)) 

TU,JME)=IT( I.JHE) 
12  CONTINUE 

IF(CKK)454,454,444 
444  CuUNT=COUNT+l. 

GO  TO  200 
454  IF(P)177,187,177 


167 
175 


185 


195 
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177  OC  +  1. 
P=P/(C0UNT+1. ) 
GO  TO  200 

187  IF(Q)178,188, 17P 

178  Q=Q/(C0UNT+1.) 
GO  TO  120 

188  GO  TO  121 
120  00  14  1=1,2 

DO  15  N=1,NN 

IF(T(I,N)-(D(N)+1. 1)240,250, 240 

240  GO  TO  15 

250  TH(I)=W(I)-SUM( I) 
IF(TH< I  1)241,251,261 

241  GO  TO  79 

251  T(I,N)=0. 

261  IF(I-2)242,252,242 

242  IA=I+1 

GO  TO  122 

252  IA=I-1 

122  IF (T( I  A, Nl- J  DIN  1  +  1. 1)243, 253,  243 

243  IF(TH(I)-(D(N)-T(IA,N) ) 1244,244,264 

244  TU,N)  =  TH< I) 
GO  TC  123 

264  T1I,N)=0(N)-T( IA.N1 
GO  TC  123 

253  IFITHI I )-D(N))245,245,?65 

245  T(  l,,\)  =  TH(I) 
GO  TC  123 

265  T(I,N)=D(N) 

123  SUM(I)=SUH(I)+T(I,N) 
15  CONTINUE 

14  CONTINUE 
121  SUMT(3)=0. 

DO  607  K=l,2 
IF(K-1)642,642,652 


104 


642  KA=K+1 

Gu  TC  6C7 
652  KA=K-1 
607  CALL  ALLT3  (T,KtH» D, NN , KA ) 

DO  16  N=1,NN 

T( 3,N)=D(N)-T(1,N)-T(2,N) 

IF(T(3,N))246,256,256 

246  GO  TO  79 

256  SUMT(3)=SUMT(3)+T(3,N) 

16  CONTINUE 
IF(W(3)-SUMT<3) ) 247, 25 7, 247 

247  GU  TO  79 

257  SUMF=0. 

DU  17  N=1,NN 

SUMC(N)=0. 

00  18  1=1,3 

F(  I,N)  =  A( I,N)»T(I,N)  +  B(I,N)«r( I,N)*<? 

SUMC(N)=SUMC(N)+F( I,N) 

18  CONTINUE 
SUMF=SUMF+SUMC(N) 

17  CONTINUE 
IF(SUMF-SMIN)248,24H,2Sg 

248  SMIN=SUMF 
JOMIN=JO 
JMIN=J 

GO  TO  258 
258  DO  19  N=1,NN 

WRITE(2,108)T(1,N),T(2,N),T(3,N) 

19  CONTINUE 
WRITE(2,109)J0MIN, JH1 N , SHIN, SUMF 

8  CONTINUE 
6  CONTINUE 

LO=LO+1 

NN=3 

IF(L0-1)921,922,921 
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922  GO  TO  900 
921  STOP 

END 

SUBROUTINE  ALLT1 (KMA, KMB, I , I  A, T,D, TH ) 

DIMENSION  TO, 4), 0(4)  ,TH(2) 
11   FURMAT(1X,11HTHIRDCHECK=,2F9.2) 
115  FORMAT! IX , 12HS ECONDCHECK=, F9 .2 ) 

Tl I,KMA)=D(KMA)-T( I A, KMA) 

W*ITE(3,115)T(I,KMA) 

IF(T(I,KMA)-TH( 1)1442,442,462 

442  T(1,KMB)=TH(I)-T(I,KMA) 
GO  TC  321 

462  T(I,KMA)=TH( I) 
T( I,KMB)=0. 

321  WRITE(3,110)T(I,KMA),T(I,KMB) 
RETURN 
ENO 

SUBROUTINE  CALHV1 LL, ZOP,A,B, I , J, HV ) 
DIMENSION  ZOP(2t8),A(3,4),B(3,4>,T(3,4) ,HV<3> 

HV(LL)=(ZOPU,J)+A(I,LL)-A(3,LL))*T(I,LL)+  i U , LL ) »T ( I , LL) »*2 
RETURN 
END 

SUBROUTINE  ALLT2 (KMA, KMB.KMC, I , I A, T,D , TH) 
DIMENSION  T(3,4),0(4) ,TH(2) 
T(  I,KMA)=D(KMA)-T(  I A, KMA) 
IFITI I,KMA)-TH( 1)1443,443,463 

443  TH(I)=TH(I)-T(I,KMA) 

T(  I,KMP)  =  0(KMB)-T(  I, KMB) 

IF  (T(I,KMI3)-TH(  I)  1444,444,464 

444  T( I,KMC)=TH(I)-T(I,KMB) 
GO  TO  322 

464  T(  I,KMB)  =  TH( I ) 
T( I,KHC)=0. 
GO  TO  322 

463  T( I,KMA)=TH( I ) 
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m,KK3)=0. 
T(  I,KKC)=0. 
322    RETU 
Ei'JD 

SUBROUTINE    ALLT3I"T , I, W, C,NN, I ' ! 
DIMENSION    T(3,4),     (?)  ,     ( ', ) 
SUMT=0. 
UJ    1    N=J  ,:.. 

1  SUMT=SUf  T  +  TU,N) 
TR=W( I >-SU«f 
IFt"TR)20>20,  L  : 

13    IF<T3-Nfmi,20,2( 

11  N=0 

2  N=N+1 
IF(T(I,f4)-D(N))12,  2,2 

12  IF((T(  it  ,N)+T(I,N)  )-Dt'i})13,2,  3 

13  IFIT(I,N))3,3,14 

3  If  (   !-A  12,14,14 

14  T  ( I ,  N )  =  T  ( I ,  '4 )  + 1 . 
TR=TS-1. 
i;:(T;}2C,?-;,2 

20    R.  TU<". 

end 

SUBROUTINE    CALZ3I  A,B,  I ,  NS,  IM,  J.H.KM 

DIMENSION    4(3,4),  V.  [3,41 

Z1=2.»3{I, IM)»B(I,JM)*  >[I,KK1*TSUM 


,JSW  ,2) 


Z2=Htifj;-:)*e{i,Kf,i»(A(  i.inj-ains.imjj 

Z3=E(t,  IM)*BU,Ktf)*U<  t,JM)-A(NS,JI  )! 
Z4=3U,IM)*DU,JM)»(A<  f,XK)-A(NS,KN]  ! 
Z5=B(I,JM)«B(I,KM)+'l(I,I,vi)«B{I,KMJ+e(I,: 

1=      (Zl  +  Z?i?J  +  /i,)/7.T 

RETURN 


}*°{I,  M) 


SUBROUTINE  CALZ2W  ,Q,I  ,IM,  JM,TSl)«,NS,Z  ) 
DIMENSION  At  3,  4)  ,:  (  i,  <,  ) 
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Z1=-2.*BU,IM)*B<I  ,JMJ*TSU  i-3<  I,  JM  )•"  (  A  (  I ,  i  ,-  1-AINS,  IK)) 

Z2=S( I,  IM)*( M I,JN)-A(NS,J  ! ) 

Z=(Zl-Z2)/(3( 1,1  ) I  (  [  ,JH) ; 

RETURN 

EN? 

SUBROUTINE  CALTKZ,  A,  ?  ,  I,IM,NS,TiI 

DIMENSION  A{3i4),3{3,4),T{3,4),ITC 

IT(I,IM)=-(Z+A(  I,  I.)- A  (OS,  IK)  )/(2.< 

Til,  IM)=-(Z+AI  I,  !N  )-,!(.<,,  I  '))/(?. 


) 
,41 

*8  < I , I M ) ! 

"  ( I ,  I  :• ) ) 


l  r-  (  t  (  ( ,  I    )  - 

in  i ,  iMi-c 

.5)1,2,2 

1   T(I,I«)=iTtl,IMJ 

CO    TC    3 

2    T I  I , IX )  =  IT  1 1 , TM  )  +  1 . 

3    RETURM 

END 

/* 
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5ULTS 


EXA«PLE  1   3  ORIGINS  MID  <  2HMAHC  POI  1" 


1  fJ=  1  .  2U 
1  N=  ?  ZU 
1  H=       3  IX! 


1  N= 

1  J= 

1  J= 

1  J  = 

I  J  = 

2 
2 


=   2 


2  J=  1 

2  J--"  2 

2  J=  3 

2  J=  A 


PI  [,.'!)=    -0.500  ZLOVl  I  ,N)  = 

»(Ii  !)=     6.000  ZLGWI  i  ,:.)  = 

MItN)=    -2.500  ZLOW(I,N)= 

-I. COO  ZLDWUtNM 


ZLV(  I,'J)  = 


Z3IJI  I,  J  )  = 
ZHUt I, J )- 
Z3UI  I,JJ  = 
7.  ill!  1,J)  = 


6 . 0  "  OC 
•0.»5000C 

•i.OOO:;C 
■2.500  )0 


1  ZUP(I,N)  = 

2  ZUP(I,N)= 

3  ZUf>U,N)  = 
<  ZUP(f,N)  = 


2.503      ZL0tf(I,N}  = 

3.500       ZLOWl 1,     )= 
■2*000       ZLOK(I,N)  = 


-0.500 

6.00  '■ 

-2.500 

-5.800 


l.r,oc 


ZL0W1  [,,-.}  = 


2  •  50r' 
3.50  : 
•3 .  OOC 

:  ■-;  _  Q  Q  Q 


Z  U( I,J)=  3.50000 
Z3UI I,J>=  2. 5C0C ( 
Z;U(I,J)=   1.00000 

zrui i, j )=  -2.C0C0C 


5.995 
-0.500 
-C.5  )5 
-1.005 

-2.505 
-5.805 
3.500 
3.495 
2.4V5 
0.995 


-2.005 

-■'  .005 
-15  i.005 
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.)>..=- 


J=    5         ZOPI  I-ltJU)=         -2.500         IQPiltJ)' 


-2.005 


r  ( l ,  h  ) 

0 

.0 

70 

.000 

000 

1 0  1 

.000 

:  t  o 

25 

.000 

0  00 

T<  ) 

,N] 

0. 

70. 

000 

35. 

000 

25. 

ooo 

J  OH 

iN= 

5 

JMIN=  5 

so.  oeoocc 

CO 

0.124979 
1.5025  ).' 

T  (  2 ,  M  ) 

50.00C 

33.00C 
2.00C 


T  (  3  ,  N  ) 

D .  " 

27.C0C 

53.000 


SKIN=  403.13  5U-F=  433.13 
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EXAMPLE    Z       3    ORIGINS 


3    DEf  *  JO    POINTS 


1 
1 

1  N= 

1  J  = 

1  J= 

1  J  = 


1       ZUPd.Nh 
?       ZUP(IiN)- 

3       ZU'J(  I  ,N)  = 


•1.500  ZLCWU.N)' 
&.  COO  ZLOWl  [  ,N)  = 
0.600       ZLOW(I,N)- 


ZBU(  [i  J)=      6.00000 

Z!  Ui  !  ,  J)  -      0.60000 

zeu(  f.J)=  -1.50000 


H-       1       ZUP(I,N)  = 

N=      2       ZOPU,H)  = 
N=      3      ZUP(I,N>= 


-1.600  ZLOWtl.Nh 
6.300  ZLOV  (I  ,N)  = 
1.600       ZLOWII,     )■ 


J=       1  ZBU(IfJ)=      6.30001 

J=       2  Z  .U(  I ,  J  )  -       1.60000 

J=       3  -ZP.U1  [,J)  =    -1.6000C 


- 1 .  r>  V  ; 
4. BOO 
0.600 


-1.600 

'..20  J 
3.80  0 


5.995 

4.795 

P.  600 

0.595 

-1.500 

-1.505 

6.300 

6.295 

!  .  5  1 0 

•  795 

-1.600 

-1.605 


JC= 


J=    2 


Z0P<I-1,J0J  = 


1.600         ZOPtIiJ)=  6.295 


f  (1,  .) 
0 . 0 


T{2,N] 


Ill 


5    .  100000 

1C 

.OO'OOi 

41.000000 

"■ 

.0 

T  (  I ,  N ) 

r(2,M) 

T{3,j*) 

0,( 

>0.< 

3O;0 30 

JO .  {     c 

0.( 

20.000 

o.a 

2     .  C    3 

JOMiN=  3    JMIN=.  P     SMIN=    452.00     SUKf=    452.00 
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The  transportation  problem  is  to  minimize  the  cost  of 
sending  a  resource  from  s  orgins,  where  the  resource  is 
located,  to  N  destinations  (demand  points)  where  the  demand 
for  this  resource  exists.   Suppose  that  there  is  one  type 
of  resource  and  that  its  total  supply  is  equal  to  the  total 
demand. 

The  aim  of  this  report  is  to  explore  the  application  of 
a  discrete  version  of  the  maximum  principle  to  the  problem 
and  to  modify  and  define  some  of  the  concepts  so  that  a 
computer  program  may  be  written  for  solving  a  problem  with 
a  large  number  of  origins  and  demand  points. 

.  The  computer  program  is  written  primarily  for  a  problem 
with  three  origins  and  four  demand  points  with  one  origin 
consisting  of  linear  cost  functions.   This  can  in  some  in- 
stances be  generalized  for  3  origins  and  N  demand  points. 

A  problem  with  three  origins  and  four  demand  points  is 
used  as  an  example  to  test  the  logic  of  the  computer  program. 


